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ABSTRACT 


This  thesis  discusses  estimating  procedures  in  Logit 
analysis.  The  classical  methods  are  outlined  and  their 
advantages  and  weaknesses  discussed.  A  new  sampling 
procedure,  the  First  Zero  sampling  procedure,  is  then 
presented.  Its  main  advantage  is  that  it  concentrates 
observations  in  the  extreme  quantile  area  of  the  response 
function,  so  it  applies  to  analyses  where  these  values  are 
the  quantities  of  interest.  The  First  Zero  distribution  and 
its  two  limiting  distributions  are  derived.  One  limit 
result  holds  in  the  general  case,  while  the  other  holds 
under  more  restrictive  limit  conditions.  How  fast  this 
convergence  takes  place  is  described  graphically.  Estimates 
of  location  and  scale  parameters  of  the  response  function 
are  derived  for  each  of  these  distributions.  The 
performance  of  these  esimates  is  compared  under  various 
sampling  situations,  and  conditions  are  determined  when 
each  of  the  estimates  offers  advantages  over  the  others. 
The  estimates  are  also  compared  to  classical  estimates  and 
it  is  found  that  in  some  cases  First  Zero  estimates  offer 
improved  efficiency  and/or  simplicity  over  the  classical 


ones  . 


t  :  ,  ,  =:  M  « 

' 

' 

■ 


ACKNOWLEDGEMENT 


One  of  the  most  important  decisions  a  student  must 
make  is  the  choice  of  a  supervisor.  In  this  I  believe 
I  was  extremely  fortunate:  Dr.  Don  McLeish  has  ful¬ 
filled  all  my  expectations  and  more.  He  provided 
insight  where  mine  was  lacking  and  suggestions  when 
the  work  had  stalled.  I  would  like  to  thank  him  for 
his  encouragement  and  patient  guidance.  I  know  I  could 
not  have  completed  this  thesis  without  his  supervision. 

I  would  also  like  to  thank  Mrs.  Sherry  Smith  for 
typing  the  mathematical  symbols  in  the  text. 


v 


. 


1 


, 


’ 


;  • 


TABLE  OF  CONTENTS 


Chapter  Page 

I.  The  Classical  Methods  . 1 

1  .  1  Introduction  . . 1 

1.2  Classical  Probit  Analysis  . 5 

1.3  Sequential  Approximation  . 12 

II.  The  First  Zero  Distribution  . 18 

2.1  Introduction  . , , . 18 

2.2  Weak  Convergence  of  e^D-1  to  Exponential  . 21 

2.3  Graphic  Comparison  of  the  First  Zero  and 

Exponential  Distributions  . r, . ...26 

2.4  Weak  Convergence  of  3  ( X  -  0 )  +  In ( ef  -  1 )  to 

Extreme  Value  . 36 

2.5  Graphic  Comparison  of  the  First  Zero  and 

Extreme  Value  Distributions  . 40 

III.  Estimating  the  Parameters  . 5] 

3.1  Introduction  . 51 

3.2  Exact  Maximum  Likelihood  Estimates  . 53 

3.3  Exponential  Estimates  . 5S 

3.4  Extreme  Value  Estimates  . 65 

IV.  Relative  Performance  of  Estimates  . 67 

4.1  Introduction  . ....6  7 

4.2  Estimates  when  3  is  known  . 70 

4.3  Estimates  when  3  is  not  known  . 73 

4.4  Performance  of  Estimates  in  Probit  Analysis  ..81 

4.5  Comparison  of  First  Zero  and  Classical 

Estimates  . 85 

V.  Conclusions  . 89 

5.1  Cost,  Loss,  and  Risk  Considerat ions  . 89 

5.2  Sequential  Approximation  Procedures  . 95 

5.3  Recommendations  . 9  7 

5.4  Further  Research  . 100 

BIBLIOGRAPHY  . 102 

APPENDIX.  APL  Programs  Used  in  Simulations  . 104 


vi 


, 


■ 


LIST  OF  FIGURES 


Figure  Description  Page 

2.3. 1- 2. 3. 9  Graphic  comparison  of  the  First  Zero 

and  Exponential  distributions  for 

various  values  of  Xj  and  A . 27-35 

2.5. 1- 2. 5. 9  Graphic  comparison  of  the  First  Zero, 

Extreme  Value,  and  Truncated  Extreme 

Value  for  various  values  of  X  ^  and  A  ...4  2-50 

LIST  OF  TABLES 

Table  Description  Page 

4.2.1  Bias  and  Variance  of  Estimates  when  3  is  Known  ..7i 

4.3.1  Bias  and  Variance  of  Estimates  when  3  is  not 

Known  and  samples  contain  10  zeros  . 7  6 

4.3.2  Bias  and  variance  of  estimates  when  3  is  not 

Known  and  samples  contain  20  zeros  . .7  7 

4.4.1  Bias  and  variance  of  Exponential  and  Extreme 

Value  estimates  in  Probit  analysis  . 83 

4.5.1  Bias  and  variance  of  estimates  of  3  and  LD95  ....87 

5.1.1  Expected  number  of  observations  per  zero  . 9l 

5.1.2  Comparison  of  WISE  of  various  estimates  of  LD95 

with  cost  restraints  . 91 


vi  1 


' 


*-1 


CHAPTER  I 


The  Classical  Methods 

1.1  Introduction .  The  problem  under  investigation 
originates  in  pharmacology,  although  it  may  be  applied  to  a 
number  of  other  situations.  The  concern  is  to  establish  the 
relation  between  a  stimulus  and  a  response.  A  common 
example  is  in  examining  the  potency  of  a  certain  drug.  We 
are  concerned  with  quantal  responses.  A  quanta!  response  is 
a  zero-one  response  (or  life-death,  cured-not  cured,  etc.). 
In  short,  any  response  that  .can  be  classified  in  such  a  way 
that  a  specific  condition  can  be  observed  as  either 
occurring  or  not  occurring  is  a  quantal  response 
experiment.  In  the  example  mentioned  above,  ie.  that  of 
examining  the  potency  of  a  certain  drug,  the  proportion  of 
the  population  that  will  respond  with  a  1  to  the  dose  z  is 
a  function  of  z.  As  z  increases,  the  proportion  of  1 
responses  increases  to  1,  and  as  z  decreases  the  proportion 
of  1  responses  decreases  to  0.  In  fact,  the  dose  itself  is 
not  usually  the  quantity  of  interest,  but  the  logarithm  of 
the  dose.  If  we  write  x=ln[z]  (either  natural  or  common 
logs  may  be  used  since  they  are  just  linear  transforms  of 
each  other),  then  we  shall  by  an  abuse  of  the  language  call 
x  the  dose.  For  a  given  dose  x,  the  proportion  of  1 
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responses  at  that  dose  is  a  function  of  x  we  wi 1 1  denbte  by 
F(x).  F ( x )  is  called  the  response  function.  It  can  be 
easily  seen  that  F  is  a  distribution  function.  Many 
functions  have  been  proposed  for  F,  but  only  two  have 
gained  any  widespread  acceptance.  If  F  has  the  form  of  a 
normal  distribution  function,  the  analysis  is  called  a 
Probit  analysis.  (Originally,  when  computations  were 
performed  by  hand,  the  term  probit  described  a  normal 
variate  shifted  to  the  right  by  five  standard  deviations. 
This  transformat  ion  avoided  the  use  of  negative  numbers. 
Although  this  is  still  done  today,  any  situation  using  the 
normal  response  function  is  termed  probit  even  though  the 
actual  "probit"  units  are  not  used.)  When  F  is  given  the 
form  of  the  Logistic  distribution,  the  analysis  is  termed 
Logit. 

At  present,  probit  analysis  is  more  commonly  used. 
This  is  mainly  due  to  Finney.  In  his  books,  "Statistical 


Method  in 

Biological 

Assay" 

and 

" Probi t  Ana lys i s " 

he 

recommends 

the  Use  of 

the 

norma  1 

distribution  as 

the 

"natural"  one.  We  quote  from  the  former :  "When  the  reason 
for  unlike  behaviour  of  similarly  treated  subjects  is 
primarily  their  intrinsic  differences  in  susceptibility, 
the  specification  in  terms  of  a  frequency  distribution  of 
individual  tolerances  is  natural,  and  ...  the  assumption  of 
a  normal  distribution  of  log  tolerances  seems  the  most 
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reasonable  procedure  in  the  absence  of  evidence  for  any 
alternative."  In  the  latter  he  does  admit:  "These  two  are 
very  similar  indeed  in  all  respects  except  for  very  small 
or  very  large  P,  and  extremely  large  experiments  would  be 
needed  to  show  one  as  a  better  fit  than  the  other.  No  one 
should  believe  that  either  formula  for  P  represents  perfect 
truth,  and  therefore  perhaps  nothing  other  than  personal 
inclination  can  decide  which  is  to  be  used."  However, 
Finney  uses  probit  analysis  exclusively  and  since  his  books 
have  become  the  authorities  in  the  area  of  bioassay,  most 
researchers  follow  suit  and  use  the  normal  as  well.  Berkson 
[  2],  however,  claims:  "In  view  of  the  wide  use  of  the 
normal  curve  to  represent  the  distribution  of  biological 
traits  and  also  because  of  direct  experimental  evidence  of 
the  normal  distribution  of  susceptibility,  it  is  to  be 
conceded  that  the  integral  of  the  normal  curve  recommends 
itself.  However,  the  logistic  function  is  very  near  to  the 
integrated  normal  curve,  it  applies  to  a  wide  range  of 
physicochemical  phenomena  and  therefore  may  have  a  better 
theoretic  basis  than  the  integrated  normal  curve.  Moreover , 
there  are  reasons  for  believing  it  to  be  easier  to  handle 
statistically."  We  concur  with  Berkson  and  in  the  methods 
we  propose,  we  assume  we  are  dealing  with  a  logit  analysis. 
For  completeness  we  include  a  brief  description  of  the 
standard  method  of  probit  analysis.  We  then  conclude  the 
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chapter  with  a  summary  of  the  Robb i ns -Monro 
which  is  the  most  commonly  used  form  of 
approximation  in  this  field. 


procedure , 
sequent i a  1 
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1.2  Classical  Probit  Analysis.  Most  of  what  is  given  in 
this  section  is  a  very  brief  summary  of  the  techniques 
outlined  in  Finney's  book  11  Probit  Analysis".  The  same 
method  could  very  easily  be  adapted  to  use  the  logistic 
distribution.  At  this  point,  however,  it  is  the  method  that 
is  of  interest  and  not  the  model.  We  will  describe  the 
procedure  used  for  estimating  drug  toxicity.  The  most 
commonly  sought  value  is  LD50,  or  Lethal  Dose  50.  This  is 
the  dose  that  would  kill  50  percent  of  the  population  under 
consideration.  (Similarly,  LD90  is  the  dose  that  would  kill 
90  percent  of  the  population.)  A  dual  quantity  which  is 
equivalent  statistically  is  the  ED50,  the  Effective  Dose 
50.  This  term  would  be  used  to  describe  a  drug's  medicinal 
potency  and  refers  to  the  dose  that  would  cure  50  percent 
of  the  population.  We  assume  that  all  trials  are 
independent,  meaning  that  no  subject  receives  more  than  one 
dose. 

If  a  subject  receives  an  administration  of  dose  X,  the 
probability  that  it  dies  (ie.  has  a  response  of  1)  is 


(1.2.1)  F ( X ) 


fX 

(  2  TT 


—  CO 


O' 


r 


exp{ - ( x- y) 2/2g2 ) }dx . 


In  this  case,  the  parameters  of  interest  are  p  and  a  (  p  is 
the  LD50).  Commonly,  the  transformat  ion  Y  =  a+PX  is  made 
making  the  parameters  a  and  3  the  parameters  of  interest 
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wherea  =-y/a  and  3=1  /a.  The  method  of  estimation  is  maximum 
likelihood.  A  typical  experiment  would  take  k  samples  (k 
would  normally  range  from  2  to  5)  of  n  subjects  each  (n  may 
range  from  5  to  500.  It  is  not  necessary  to  have  the  same 
number  in  each  sample  but  it  does  simplify  calculations) 
and  administer  each  subject  in  the  i th  sample  with  dose  X.. 
For  example,  if  k=3,  typical  values  for  X-^X^,  and  X^  would 
be  chosen  hopefully  near  LD20,  LD50,  and  LD80  respectively. 
The  number  r-  of  1  responses  in  the  i th  sample  is  observed. 
r_j  has  a  binomial  (n,  P  . )  distribution  where  P^=F(X^).  The 
joint  density  of  the  r. ' s  is 


i  ( 1  -  P  i ) 


n  -  r 


i 


leading  to  a  log  likelihood  equation  of 


InL  = 


=  ■  V  {r.ln[Pj  +  (n-r.  )  ln[  1-P.  ] }  . 
i=i  1  1  1  1 


(Note  that  the  n' s  in  the  above  will  be  subscripted  by  i  if 
the  sample  size  varies  from  sample  to  sample.) 

To  maximize  InL.  we  take  the  partial  derivative  of  InL  with 
respect  to  a  and  3,  set  these  equal  to  0,  and  solve.  Ie. 


9 1  nL 

9  GL 


k 

=0=1  {r  /? 
i  =  1  i  i 


(n-r  ) / ( 1-P  )}  3Pi 
i  i 


9a 
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( 1.2.2) 


91  nL 

93 


0  = 


i  Vpi 

i  =  1  1  1 


(n-r.  )/(1-P.  )} 


iLLi 

93  ' 


These  cannot  be  solved  explicitly,  but  a  two  parameter 

iteration  procedure  can  be  used  to  estimate  a  and  3. 

Suppose  that  a^and  3  are  approximations  of  a  and  3  . 

Second  approximations  to  a  and  3  will  be  of  the  form 

a,  +  Aa  and  3  +  A3  ■  Using  the  Tay lor-Macl aur in  expansion 
1*1  - 

of  InL,  these  second  approximations  must  satisfy  both  of 


91  nL 

3  2  1  rt  L 

+  A3 

9  2  1  n  L 

9a 

+  A  a 

da1 

9  a  9  3 

=  0 

9 1  n  L 

+  A3- 

9  2  1  n  L 

+  Aa 

92  1  nL 

=  0 

93 

93  2 

9a  9  3 

These  two  linear  equations  can  be  solved  for  Aa  and  A3  . 
Finney  recommends  simplifying  the  second  order  partial 
derivatives  in  (1.2.3)  by  substituting  P.  =r.  /n  in  these 

ii 

derivatives  after  differentiation,  giving  expected  values 
rather  than  empirical.  The  next  step  should  now  be  clear. 

In  general,  it  is  best  to  have  the  iteration  continue  until 

/\  ^ 

convergence  occurs.  The  covariance  matrix  of  a  and  3  is 
asymptotic  to 


I 


f  9 2 1 nL 


9a93 


V 


v- 


3  2  1  n  L 

9a  93 


-E 


9  y  "  ' 


and  confidence  bounds  can  be  set  for  a  and  b  (the  estimates 
of  a  and  3  respectively).  The  above  applies  for  any  form 
that  the  Rj  '  s  may  take.  Specifically,  in  probit  analysis, 
P  .j  has  the  form  (1.2.1),  and  if  we  wr  i  te  cf>  ( x )  and  $  ( x )  as 
the  standard  normal  density  and  distribution  functions 
respectively,  we  have  (recalling  the  transformation 


Yi  =  a  +  3  X:  ) 


Thus,  (1.2.3)  in  probit  analysis  becomes 


k  k  k 


(  1.2.4) 


k 


k 


k 

l  nw.  X.  (  r.  /n-P.  )  /  <j>(  a+bX.  ) 
i  =  l  1  1  1 


where  w,*  =  <f>2  ( a+bX  • )  /  [  P  _•  ( 1  -  P  • )  ]  . 


If  a  and  b  are  the  present  estimates,  then  by  solving  these 
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equations  for  A  a  and  Ab  the  next  estimates  a'  and  b'  may  be 
calculated  by 

a'  =  a  +  A  a  and  b'  =  b  +  Ab. 

The  procedure  continues  in  this  fashion  until  convergence 
occurs . 

Usually  a  Chi -squared  test  is  done  to  check  the 
goodness  of  fit  of  the  observed  and  estimated  data.  There 
is  a  reasonably  high  risk  of  having  a  small  expected  value 
in  some  of  the  extreme  samples  which  could  cause  problems 
in  the  validity  of  the  test. 

The  above  analysis  has  proved  itself  over  the  several 
years  since  Finney  first  presented  it  in  1947  (in  fact  the 
method  had  been  available  as  early  as  1935.  See  [4],  [5], 
and  [6]).  One  advantage  it  has  is  that  samples  where  the 
proportion  of • responses  is  either  0  or  1  do  not  affect  the 
procedure.  (It  is  possible  to  have  unreliable  results  if  a 
large  number  of  the  samples  fall  into  this  category,  but 
this  should  not  happen  frequently  in  actual  practice.)  It 
is  also  not  necessary  to  have  equal  sample  sizes,  and  there 
are  no  theoretical  problems  arising  from  having  a  sample 
size  of  one.  This  is  important  to  understand  since  the 
method  has  been  erroneously  criticized  on  both  of  these 
counts  [  8]  .  A  possible  problem  which  might  be  encountered 
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is  the  sensitivity  of  the  convergence  scheme  to  the  choice 
of  starting  estimates  for  the  parameters.  Finney  devotes  a 
large  section  of  his  book  to  getting  initial  estimates 
graphically.  This  was  practical  when  computations  were  done 
by  hand  and  accurate  initial  values  shortened  the 
procedure.  But  now  with  the  widespread  availability  of 
computers,  it  is  cumbersome  to  have  to  rely  on  'eyeball' 
procedures  to  get  the  initial  values.  The  method  does  have 
other  drawbacks  in  certain  situations,  which  we  will 
mention  later  when  we  compare  this  method  to  both  the 
Robb i ns -Monro  procedure  and  the  method  we  propose. 

We  have  presented  here  the  classical  probit  analysis. 
The  theory  is  the  same  for  logit  analysis  except  for  two 
major  differences.  The  first  is  the  further  simplification 
possible  because  of  the  form  of  the  logistic  distribution, 
which  is 


F  (  x)  =  1/(  l+e'x)  . 


Making  the  transformation  Y  =  6( X -  e)  gives 


Pi  =F(  Ti’  •  ^j=-epi(1-pi>'  and  !!L=(X.-e)P.  (1-P.). 


By  defining  w.  =  P.  (1-P..)  the  cor  respondi  ng  form  of  (1.2.3) 
in  logit  analysis  is 


' 

V, 


■ 


k  k  k 

A6  l  6  w  .  +  Ag  l  (X.  -  e)w.  =  -  )  (r.  /n-P.  ) 

i  =  l  1  i  =  1  1  1  i  =  1  1  1 


(1.2.6) 


k  k  k 

-  A0  I  (X.-0)w.+A$£  (X.-b)2w.  =  }'  (r./n-P.  )  ( X.  - 

i  =  1  1  •  1  i  =  1  1  1  i  =  1  1  1  1 


The  second  advantage  is  that  from  (1.2.6) 
apparent  that 


1  1 


it  is 


are  sufficient  statistics  for  e  and  3  respectively.  No 
sufficient  statistics  are  forthcoming  in  probit  analysis. 
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1.3  Sequential  Approximation.  This  method  of  approximation 
gets  its  name  from  the  type  of  procedure  it  requires.  In 
the  case  of  establishing  drug  potency,  it  would  be 
necessary  to  make  an  administration  at  a  certain  level  and 
then  observe  the  response  before  the  next  administration  is 
made.  This  is  because  the  next  dose  is  a  function  of  the 
previous  dose  and  the  response.  In  general  terms,  we  wish 
to  estimate  the  root  of  a  function  R(x)  when  the 
observations  are  subject  to  random  error.  Hence,  what  is 
observed  is  R(x)+  c  where  e  is  a  random  variable  whose 
distribution  may  depend  on  x.  In  probit  or  logit  analysis, 
R(x)  has  the  form  F(x)-0.5  (in  the  case  of  estimating  the 
LD50)  where  F(x)  is  the  normal  or  logistic  distribution 
function  respectively.  Robbins  and  Monro  [  14]  proposed  a 
sequential  method  that  would  define  X}  arbitrarily  and  then 
define 

(1.3.1)  Xn+1  =  Xn  -  an(R(Xn)  +£  ) 

where  an  is  a  positive  decreasing  sequence  of  constants 
such  that 

00  00 

l  a  diverges  and  J  an 2  converges. 
n=  1  n  n=l 

They  showed  that  when  R  increases  at  most  linearly  and  the 


> 
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variance  of  the  e's  is  uniformly  bounded,  then  X  converqes 

n 

to  the  root  of  R(x)  in  probability.  Blum  [  7]  later  showed 

that  in  fact  converges  to  the  root  almost  surely.  Since 

then  there  have  been  numerous  papers  that  have  weakened  the 

requirements  on  R  and  have  dealt  with  various  applications 

of  the  method  to  other  situations,  for  example  Kiefer  and 

Wolfowitz  [  13]  have  altered  the  method  to  allow  the  maximum 

or  minimum  of  a  function  to  be  approximated  sequentially. 

Others  have  considered  the  problem  of  R  having  multiple 

roots.  A  considerable  amount  of  effort  has  gone  into 

finding  the  optimum  values  for  the  sequence  {an}.  For 

example,  it  has  been  shown  that  among  all  sequences  of  the 

form  c/na,  the  optimal  form  is  a  =c/n. 

n 

This  procedure  when  compared  to  classical  analysis  has 
advantages  and  disadvantages.  Its  primary  disadvantage  is 
its  sequential  nature.  When  testing  the  potency  of 
insecticides  it  is  highly  undesirable  to  fumigate 
individual  cockroaches  with  individual  doses  and  then 
observe  the  life  or  death  response  before  the  next  test  is 
performed.  Cockroaches  are  cheap  and  the  cost  of  testing 
high  so  in  this  situation  sequential  analysis  is  not  cost 
effective.  Similarly,  in  testing  the  carcinogenic  effect  of 
a  new  drug  on  mice,  it  is  again  highly  undesirable  to  wait 
the  substantial  period  of  time  required  for  the  effects  to 
become  noticeable  before  the  next  administration  is  made. 


' 

■ 
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In  this  case  the  method  is  not  time  effective. 

The  method,  however,  does  have  several  significant 
advantages.  Chief  of  these,  perhaps,  is  its  simplicity. 
Finney  when  describing  the  use  of  probit  analysis  without  a 
computer  required  several  pages  of  explanation  as  well  as 
some  extensive  examples.  The  sequential  procedure,  on  the 
other  hand,  could  be  done  by  hand  or  small  pocket 
calculator  right  in  the  laboratory  by  someone  with  a 
minimum  of  theoretical  or  computational  skill.  Also,  there 
are  many  experimental  situations  involving  expensive 
animals  (for  example  Rhesus  monkeys)  and  fast  acting  drugs 
that  make  the  method  very  cost  and  time  efficient.  Another 
situation  where  this  method  is  cost  and  time  efficient  is 
in  the  testing  of  plastic  pipe  for  breaking  strength.  Here 
a  given  weight  (the  dose)  is  dropped  from  a  fixed  height  on 
a  sample  piece  of  pipe,  and  the  response  of  breaking  or  not 
breaking  is  recorded.  The  very  nature  of  the  testing 
procedure  lends  itself  very  nicely  to  sequential 
approximation.  Another  of  this  method's  definite  advantages 
is  that  it  is  nonparamet r i c .  The  ability  to  estimate  and 
predict  various  LD  values  in  probit  or  logit  analysis 
relies  heavily  on  the  underlying  distribution  that  is 
chosen.  In  fact  the  estimates  themselves  depend  upon  the 
assumption  of  the  model  for  their  validity.  The 
Robbi ns-Monro  procedure,  however,  will  converge  to  the  LD50 


* 
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regardless  of  the  underlying  distribution  function.  It  also 
has  the  added  advantage  of  placing  most  of  the  observations 
near  the  root  of  interest.  This  prevents  the  rather 
wasteful  occurrence  of  taking  a  batch  of  observations 
nowhere  near  the  root  and  then  finding  either  0  or  100 
percent  response  in  that  group.  As  well,  the  experimenter 
has  the  choice  of  making  the  experiment  consist  of  a  fixed 
number  of  observations  or  letting  it  continue  until  the 
desired  level  of  accuracy  is  achieved. 

For  the  most  part  we  have  referred  to  the 
Robb i ns -Monro  procedure  as  it  is  used  to  estimate  the  LD50. 
One  reason  for  this  is  that  this  quantity  is  of  great 
interest  among  researchers  today  and  is  used  extensively 
for  comparing  drugs.  Another  reason,  however,  is  that  the 
method,  although  theoretically  asymptotically  unbiased, 
demonstrates  considerable  bias  when  used  to  approximate 
noncentral  quantities  (for  example  LD75  or  LD95)  in  finite 
samples.  Wetherill  has  considered  this  problem  [16]  and 
concludes  that  "Routine  1 ' ( the  Robbi ns-Monro  procedure)  is 
very  efficient  for  estimation  of  L.5  (ie.  LD50),  both  as  a 
method  of  placing  observations  and  as  a  method  of 
estimation.  It  is  very  robust  to  errors  in  the  starting 
value  of  the  sequence  and  also  to  the  value  of  the  constant 
c  (recall  an=c/n).  Actual  small  sample  variances  follow 
closely  approximations  given  by  asymptotic  formulae  based 


' 

<  : 


16 


on  a  simple  model.  However,  Routine  1  is  unsuitable  for 
estimating  even  moderately  extreme  Lp,  such  as  L.75,  being 
subject  to  a  heavy  bias  and  yielding  large  variances." 
Finney's  method  must  also  be  altered  for  the  case  when 
extreme  quantiles  are  to  be  estimated,  and  he  recommends 
using  a  sequential  method  (where  practicable)  proposed 
originally  by  Bartlett  [  1  ].  This  involves  sampling  at  a 
fixed  dose  until  the  desired  number  of  responses  is 
observed,  thus  yielding  a  negative  binomial  rather  than 
binomial  distribution  for  the  number  of  zeros  at  each  dose. 
Tsutakawa  [15]  has  also  proposed  a  method  for  determining 
what  the  optimal  choices  of  dose  would  be.  Wetherill  [17] 
has  suggested  a  Robbins  -  Monro  -  like  procedure  which 
would  have  smaller  bias  of  extreme  quantile  estimates.  His 
UDTR  (Up  and  Down  rule  on  a  Transformed  Response  curve)  is 
a  variation  of  the  up  and  down  procedure  first  put  forth  by 
Dixon  and  Mood  [  9  ]. 

The  estimation  procedure  we  propose  originated  from 
studying  this  particular  problem.  Our  aim  was  to  find  a 
procedure  which  would  efficiently  estimate  extreme 
quantiles  such  as  LD75,  LD95,  or  LD99.  Such  procedures  are 
becoming  important  in  areas  such  as  environmental  control 
and  establishing  safe  levels  for  toxic  drugs.  Vie  feel  there 
is  a  need  for  a  method  that  provides  high  accuracy  with 
fewer  0  responses  than  the  classical  methods.  In 
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particular,  our  methods  will  provide  the  same  accuracy  of 
estimates  with  fewer  deaths  in  situations  such  as  testing 
the  polluting  effects  of  certain  chemicals  on  wildlife. 


CHAPTER  II 


The  First  Zero  Distribution 

2.1  Introduction .  We  decided  to  call  the  distribution  which 
we  concentrated  on  the  First  Zero  Distribution.  We  assumed 
we  were  dealing  with  the  logistic  model  (or  Logit)  which 
has  the  distribution  function 

F ( x )  =  1  / ( 1  +  e'B(x~6)  ). 

3  is  the  scale  parameter  (  3  is  inversely  proportional  to 
the  standard  deviation)  and  0  is  the  location  parameter  (e 
is  the  mean  of  the  distribution).  In  the  model,  subjects 
are  administered  various  doses  of  a  drug.  The  subjects  are 
observed  for  either  of  two  possible  reactions  -  response 
vs.  non-response,  death  vs.  survival,  etc.  -  that  are 
classified  as  either  0  or  1  .  We  will  denote  doses  by  x  and 
responses  by  Y.  Hence,  at  any  given  dose  x,  the  probability 
of  having  a  response  Y  =  1  is  F(x).  We  are  interested  in 
determining  relatively  large  dosage  levels,  for  example 
LD95.  The  test  procedure  goes  as  follows:  Choose  an 
arbitrary  starting  point  X  (hopefully  chosen  near  LD95) 
and  an  increment  value  a  .  The  first  subject  is  administered 
dose  X  and  the  reaction  observed.  If  the  response  Y^  is  0 
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then  N  =  1  and  we  estimate  LD95  accordingly.  (N  will  be  the 
random  variable  correspondi ng  to  the  number  of  subjects 
treated  before  a  0  is  observed.)  If  the  response  Y  is  1, 
then  we  set  =  X  -  A  and  the  next  subject  is 
administered  dose  X^.  If  the  response  is  0,  then  we  set 
N  =  2  and  the  estimate  of  LD95  is  calculated.  If  the 
response  Y?  is  1,  then  we  set  X  =  X  -  a  and  repeat.  The 
testing  continues  in  this  fashion  until  a  0  is  observed 
(whence  the  name  First  Zero  Distribution).  Since  we  are 


stepping  to 
The  estimate 
def i ned : 


the  left,  a  0  must  occur  with  probability  1. 

is  based  on  X  ,  X  ,  and  A.  Hence,  N  may  be 

I  N 


N  =  inf  { j  :  Y.  =0}. 

J  J 

Equivalently,  we  may  base  the  estimate  on  X  ,  a,  and  D 
where  D  is  the  random  variable  defined  by 
D  =  =  a(N-1). 

The  resulting  distribution  can  be  explicitly  written. 

Assuming  the  subjects'  responses  to  be  independent,  the 

probability  of  a  response  of  Y  =1  from  the  i'th 

i 

admi ni s traion  is 


p(Yi  =1  )  =  F(Xi  )  =  1/(1  +  e“ S (Xi  -  e )  ) 


where  X-j  is  the  dose  administered  to  that  subject.  Hence, 
the  probability  that  none  of  the  first  n  responses  are  0  is 


O' f 
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P(X  <X  )  =  n  1/(1+  e  B(xre)  ) . 

N  n  i  =  j 

Although  we  have  the  explicit  First  Zero  Distribution, 
it  does  not  lend  itself  to  easy  manipulations,  especially 
in  the  area  of  parameter  estimation.  We  had  hoped  to  find 
simpler  distributions  which  would  approximate  the  First 
Zero  Distribution  and  thereby  facilitate  the  estimation 
procedure.  In  this  we  were  successful:  as  goes  to 
infinity  and  A  goes  to  zero  a  function  of  D  has  a  limiting 
Exponential  distribution. 


' 
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2.2  Weak  Convergence  of  e  1  to  Exponential.  To  help 
simplify  notation  we  introduce  the  quantity 

(2.2.1)  m  =  e'e(xr0V(eSA-1  )  . 

Theorem  2.2.1:  Let  the  constant  z>0  be  given.  Then 


P(m(epD-1 )  <  z)  -  1 -e" z 


as 


Proof : 


X  ->oo  and  a 


Define  d  and  n  by 


0. 


(2.2.2)  d  =  ( ln( 1+z/m) )/3 


and 

(2.2.3)  n  =  1  +  l d/ A 1 

=  1  +  greatest  integer  part  of  d/A . 


Note  that 


(2.2.4)  0  <  nA-d  <  A 
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and  also  that 

(2.2.5)  z  =  m(e^^-l). 

Now  P(D>d)  =  P(Y-  =1  ,Y9  =1 . Y  =1  ) 

1  2  n 

=  n  1/(  1+e  “eM 
i  =  1 

and  hence 

(2.2.6)  -In  P(D>d)  =  J  1  n  (  1  +e  ‘  B  ( xi  ~e  ))  . 

i  - 1 

Using  the  inequality  a-a2  <  ln(1+a)  ^  a  for  all 

termwise  on  (2.2.6)  gives 

(2.2.7)  l  e-p(Xi -0)  -  I  e-2B<xi-6) 

i  =  1  i  =  l 

<  -In  P ( D>d ) 

<  ie"e(X  i-.e) 

i  =  1 

which  factors  into 

(2.2.8)  e-S(Xr9’nj1eeAi  -  e-ze(X1-0)nj1e26A1 

i =0  i =0 


a>0 
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<  -In  P ( D>d ) 

<  e-f3(xr6)ni1e6Ai. 

1=0 

Summing  the  geometric  series  gives 

(2.2.9)  e'e(xr0)(eenA-1  )/(e6A-1  ) 

-  e-20^r6)(e26nA-1)/(e2eA-i) 

<  -In  P ( D >d ) 

<  e-B(Xj-6)(e3nA_'|)/(e0A-i) 
which  is  equivalent  to 

(2.2.10)  m(e<3nA-1)  -  m2  ( e2  0nA- 1  )  ( e0A- 1  )  /  (e3A+ 1  ) 

<  -In  P ( D>d ) 

S  m  ( e0  n  A  -  1 )  . 

We  wish  to  show  that  the  second  term  in  the  left  hand  side 
of  (2.2.10)  goes  to  zero  as  Xj  tends  to  infinity  and  A 
tends  to  0;  ie.  we  wish  to  show  that 
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(2.2.11)  m2(e2f3nA-1  )  (eeA-1)/(e6A+1  )  -  0. 


To  see  this,  note  that 

0  <  m2(e2(3nA-1  )  (eeA-1  )/(eea+1  ) 

^  m2  ( e 2^nA- 1 )  (e  ^A- 1 ) 

<  m2 (e2^( d+A )- 1 ) (e BA- 1 )  [from  (2.2.4)] 

^  m2 ( e28A( i+z/m) 2 - 1 ) ( e ®A- 1 )  [from  (2.2.2)] 

£  m2(e2B&-l  +  2ze2BA/m  +  z2e  23*/m2 )  (e  1 ) 

<  (eeA+1)e-2e(xr0)+2ze-B(xr6-2A)+z2(eBA-1)e2B& 

[from  (2.2.1)] 

The  first  two  terms  go  to  0  as  tends  to  infinity.  The 
last  term  goes  to  0  as  A  tends  to  0. 

Next  we  wish  to  show  that 

(2.2.12)  (eBnA-  1 )  /  (eed- 1 ) ->  1 

« 

This  follows  from 

1  <  (e8nA-1)/(eed-1) 


[from  (2.2.4)  ] 


•- 
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[from  (2.2.4) ] 

[ from  (2.2.2)] 

The  left  term  tends  to  0  as  X  tends  to  infinity  and  the 
right  term  tends  to  1  as  a  tends  to  0.  Now  noting  that 

P ( D>d )  =  P(m(e3D-1 )  >  m(eed-1 ) ) 

=  P(m(e^D-1)  >  z)  [from  (2.2.5)] 

we  have,  by  applying  (2.2.11)  and  (2.2.12)  to  (2.2.10), 
that 


<  (e3(d+A)-i )/(e3d-i ) 


=  (e3A ( 1+z/m) -  1 ) / (z/m) 


=  m(e3A- 1 ) /z  +  e3A 


=  e 


"3(xre}z  + 


3A 


-  In  P (m(e3D- 1 )  >  z )  ^  z 
which  is  the  result  we  wanted. 
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2 . 3  Graphic  Comparison  of  the  First  Zero  and  Exponential 

Pi str ibut ions .  The  above  theorem  guarantees  convergence  of 
m(e^D-1)  to  the  Exponential  distribution  in  limiting 
situations.  For  practical  purposes  it  is  important  to  Know 
how  fast  this  convergence  takes  place.  The  following  graphs 
indicate  the  kind  of  approximation  that  can  be  expected  for 
6=  0,  3  =  1,  and  various  values  of  Xj  and  A  .  As  can  be 
seen,  the  approximation  is  reasonably  good  for  even  small 
values  of  X j  (ie.  values  as  central  as  L D 7 5 )  and  large 
values  of  A  (ie.  values  as  large  as  .5/3).  As  expected  the 
fit  improves  as  either  either  X  gets  larger  or  a  gets 
smaller.  When  both  occur,  a  very  good  approximation  is 
obtained  (eg.  X  =  LD99  and  a  =  . 02/ 3 ) -  The  lack  of  fit  for 
either  X^  central  or  a  large  seems  to  be  primarily  due  to 
the  situation  of  approximating  a  discrete  distribution  by  a 


continuous  one. 
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FIGURE  2.3.3 
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FIGURE  2.3,4 
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FIGURE  2,3.5 
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FIGURE  2.3.6 
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FIGURE  2.3,7 
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FIGURE  2.3.8 
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FIGURE  2.3,9 
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2.4  Weak  Convergence  of  B  (  -  o)  + 1  n  ( e  ^A-  1  )  to  Extreme  Value. 

The  results  of  section  2.2  hold  regardless  of  the  relative 
rates  of  convergence  of  X  to  infinity  and  A  to  zero.  If  in 
fact  this  relative  convergence  takes  place  under  certain 
conditions,  another  limiting  distribution  can  be  shown  to 
apply.  Recall  that 

(2.4.1)  D  =  X2  -  XN 
and 

(2.4.2)  m  =  e~G  (^i"G  V  (eGA- 1  )  . 


Corol  1  ary  2.4.1:  Let  the  constant  z  be  given.  Then  as  X  ]_ 
tends  to  infinity  and  A  tends  to  0  in  such  a  way  that  m 
tends  to  0  we  have  that 


P(  3(X^  "  e)  +  ln(eGA-  1  )  <  z)  ->  exp(-e"z). 

Proof :  Putting  x=e"z  we  have  that 

1-e"x  =  lim  P(m(e^D-1)  ^  x)  [from  Theorem  2.2.1] 
=  lim  P(me^D  ^  x)  [since  m 0  ] 
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=  lim  P(me^Xl  XN^  <  x)  [from  (2.4.1)] 

=  lim  P(e"3(XN~0)/(e3A-1)  <  x)  [from  (2.4.2)] 

=  lim  P(  -  3(X^|- 6)  -  ln(e^A- 1  )  ^  -z)  . 

Hence  P  (  3  (  X^- 0 )  + 1  n  ( e  ^A-  1  )  <  z)  ->  exp(-e“z). 

exp(-e“z)  is  the  standard  form  of  the  Extreme  Value 
distribution  (sometimes  called  Gumbel's  Extreme  Value 
distribution  or  the  Type  I  Extreme  Value  distribution). 

To  understand  the  implications  of  this  result  we 
should  look  more  closely  at  the  condition  that  m-*  0.  From 
(2.4.2)  and  using  the  fact  that  ( e 3A-  1  )  /  3A+  1  asA  ->  0  we 

see  that 

m^  0  iff  Ae3Xl  iff  3  X^  +  In (  a) 

These  equivalent  conditions  indicate  more  clearly  the 
relative  sizes  of  X  and  a  .  In  general  terms  the 
requirement  may  be  expressed  by  stating  that  a  may  not 
converge  to  0  too  quickly. 

One  further  consideration  should  be  made.  In  practice 
Theorem  2.2.1  will  be  interpreted  to  mean  that  since 


I 
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lim  P  ( m  ( e  ^  D-  1  )  <  z )  =  1  -  e  ~  z 

we  may  conclude,  for  finite  values  of  X  and  nonzero  values 
of  A,  that 

(2.4.3)  P(m(e6D-1)  <  z)  £  1-e'z. 

Using  (2.4.1)  we  can  restate  (2.4.3)  in  the  equivalent  form 

(2.4.4)  P(X^<z)  =  exp- {exp- $(z- e) -exp- $(X  - e) }/ (e ^A- 1 ) 
or  more  simply 

(2.4.5)  P(X^<z)  =  exp- {exp- 3(z- 1 )  -  exp-ptX  -t)} 
where 

(2.4.6)  t  =  0  -  (ln(e3&-1))/e. 

This  is  the  form  of  a  Truncated  Extreme  Value  random 
variable  which  is  truncated  above  at  Xj .  That  is  (2.4.5) 
represents 

P  (  X  <z  |  XN<X  -,) 

N  1 


(2.4.7) 
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for  an  Extreme  Value  random  variable  with  location 
parameter  t  and  scale  parameter  3  . 
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2 . 5  Graphic  Comparison  of  the  First  Zero  and  Extreme  Value 

Pi str ibut ions .  The  above  theorem  guarantees  convergence  of 
3 ( Xn “ e )  +  ln(e^A-1)  to  Extreme  Value  in  certain  limiting 
situations.  Again,  we  have  graphed  the  two  distributions  on 
the  following  pages  to  indicate  the  Kind  of  fit  that  can  be 
expected  for  finite  cases.  As  for  the  graphs  comparing  the 
Exponential  and  First  Zero  distributions,  we  have  assumed  0 
and  3  to  be  0  and  1  respectively.  As  can  be  seen,  the 
approximation  is  closer  for  large  values  of  and  large 
values  of  A  (which  yield  a  small  value  of  m) .  In  the  best 
fitting  graph,  X^  is  equal  to  the  LD99  and  A  is  equal  to 
.5.  Thus  the  corresponding  value  of  m  is  0.0156.  Compare 
this  to  the  fit  for  the  cases  X^=LD99  and  A  =.02  (giving 
m=.5000),  X  =LD75  and  A  = .  5  (giving  m=.5138),  and  X^LDYS 
and  A  =.02  (giving  m=16.50).  It  becomes  quite  evident  that 
the  fit  deteriorates  rapidly  as  m  increases. 

In  order  to  compare  this  approximation  to  that 
obtained  from  the  exponential  form  we  have  on  the  same 
graph  included  the  Truncated  Extreme  Value  approximation  of 
(2.4.5).  As  expected,  in  all  cases  the  Truncated 
approximation  is  better  than  the  nontruncated  one.  Only  in 
two  cases,  namely  Xi=LD95  and  LD99  and  a  =.5  is  the 
nontruncated  fit  competetive.  We  would  therefore  expect  the 
Extreme  Value  approximation  to  have  limited  applications. 
We  shall  see,  however,  in  certain  situations  this 
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approximation  does  prove  surprisingly  useful. 
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FIGURE  2.5.1 
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FIGURE  2.5.2 
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FIGURE  2.5.3 
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FIGURE  2.5.4 
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FIGURE  2.5.5 
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FIGURE  2.5.6 
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FIGURE  2.5.7 
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FIGURE  2.5.8 
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FIGURE  2,5.9 
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CHAPTER  III 


Estimating  the  Parameters 

3.1  Introduct ion .  The  main  reason  for  approximating  a 
distribution  is  to  facilitate  the  process  of  estimating  the 
parameters.  We  will  be  concerned  with  three  parameters:  3 
(or  perhaps  1/3  ),  0  ,  and  r  where  r  is  a  p-fractile  other 

than  9.  For  example,  if  p=.75  then  r=LD75.  If  p=.5,  then 
r=LD50=  0.  Since  the  logistic  distribution  is  a  two 
parameter  distribution,  any  two  of  these  parameters 
determine  the  third.  Usually  of  interest  are  the  scale 
parameter  and  one  location  parameter.  We  will  also  consider 
the  simpler  problem  of  estimating  e  and  r  when  3  is  Known. 
There  are  three  sources  of  estimates:  the  Exact  First  Zero 
Distribution,  the  Exponential  approximation,  and  the 
Extreme  Value  approximation.  In  addition,  for  each 

distribution,  there  are  various  types  of  estimates.  These 
include  maximum  likelihood  (ML)  estimates,  moment 

estimates,  and  uniformly  minimum  variance  unbiased  (UMVU) 
estimates.  Thus  if  we  talk  about  the  Exponential  UMVUE ,  we 
are  referring  to  the  UMVU  estimate  obtained  from  the 
Exponential  approximation.  Note  that  this  is  not  a  "true" 
UMVUE.  It  is  UMVU  for  the  Exponential  distribution  but 
undoubtedly  bias  will  be  introduced  from  the  finite 
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approximation.  By  the  same  token,  the  Exponential  NILE  will 
be  asymptotically  unbiased  with  respect  to  estimating  its 
own  parameters  but  may  not  have  this  property  with  respect 
to  the  First  Zero  parameters. 

Each  of  the  next  three  sections  in  this  chapter  is 
devoted  to  deriving  estimates  from  each  of  the 
distributions:  section  3.2  covers  estimates  from-  the  Exact 
First  Zero  distribution,  section  3.3  discusses  Exponential 
estimates,  and  section  3.4  describes  Extreme  Value 
estimates.  Their  relative  performance  will  be  evaluated  in 
the  next  chapter . 


■ 
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3.2  Exact  Maximum  Likelihood  Estimates.  Given  the  explicit 
First  Zero  distribution,  it  is  possible  to  obtain  the  Exact 
Maximum  Likelihood  Estimates  for  e  and  3  (or  for  r  and  3  ). 
The  distribution  function  for  each  0  observed  is 


P(XN  =x 


n 


)  =  P(  Y1=1  f  Y2=1 


Y  =1 
n-  1 


Y  =0) 
n 


=  F(X)F(XJ...F(X  )  [  1  -  F  ( X  )] 

1  2  n-  1  n 

=  {  n  F ( X . ) }exp- 3 (X  -  0 ) 

1=1  1  n 

where  F(x)  is  the  logistic  distribution  function 
F ( x )  =  1 / [ l+exp-3 ( x- e ) ] . 


For  k  such  iterations  with  starting  points 


X-i  ,  X 1  Xi  , 

A1  Ak 

increment  sizes  A^,  •••  >  A^, 

and  first  zeros  at  X  ,  X  .  X 

n  l  n2  nk 

the  likelihood  function  is 
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L  ■  p|x«,  ■  v  \  ■  % . \ 


►  ”j 


=  {.n  n  F(x  )  }exp- 3  y 
O-i  1  =  1  1.  .1  =  1 


(X  -  0) 
J-i  J 


This  leads  to  the  log  likelihood  function 


X  ) 

nk 


1  n  (  L  ) 


k  "j 

=  I  I  ln( F ( X  .  )) 
j=l  i=l  'j. 


Taking  partial  derivatives  with  respect  to  0  and  3  gives 


31 n ( L  )  .  v  . 


ae 


=  -l  p  r  F<X  lexp-  g( x  -  0) 

j=i  i  =  i  .i  Vi 


k  e 


and 


31  n  (  L )  _ 

03 


n 


l  £JX.  F  (  X  .  )exp-3(X. 

j  - 1  i  - 1  Vi  nj  1 


-  6  ) 


k 

I 

j*l 


(X 


n 


0) 


Setting  these  equal  to  0  and  reducing  leads  to  the 
equations 


k  nj  k 

(3.2.1)  l  l  F ( X .  )  '  l  (n  -  1)  =  0 

j  =  l  i  =  1  nj  j=l  J 


and 


n 


(3.2.2)  l  IJX..  F  (X  ) 

•  -  •  1  I  •  I  * 

J=1  l  =1  J  J 


k  nj 

"I  l  xi 

j=l  i=l  J 


0. 
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These  cannot  be  solved  explicitly  but  are  manageable 
numerically.  The  standard  procedure  is  to  do  a  two 
parameter  Newton-Raphson  recursion.  If  3  is  Known,  then  a 
one  parameter  recursion  may  be  used  on  (3.2.1)  to  find  the 
value  of  0  . 

In  practice,  the  values  of  and/or  A  will  be  changed 
during  the  course  of  the  experiment.  However,  for  the  sake 
of  the  comparison  we  will  be  doing  in  the  next  chapter,  we 
also  consider  the  case  with  the  further  restrictions  that 
all  the  starting  points  and  all  the  increment  sizes  are 
equal.  Under  these  conditions  the  process  is  simplified 
noticeably  and  by  defining 


n  =  max  {n  ,  n^ ,  .  .  ,  , 


we  can  then  define  the  values  for  1  <  i  <  n 


w  =  the  number  of  n  /  s  that  are  ^  i 
i  J 


where 


I  In .  >  i 1 
J 


C  1  if  n  .  >  i 
J 


0  if  n.  <  i  . 
0 


1 


-  .. 
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This  reduces  (3.2.1)  and  (3.2.2)  to 


n  k 

(3.2.3)  I  w.  F ( X .  )  -l  (n.  -  1 )  =  0 

i=i  1  j=i  J 


and 


k  n 

(3.2.4)  l  Xn  -  l  w.  X.  (  1  -  F  ( X .  )  )  =  0. 

...  n  .  . L ,  ii  l 

J=1  J  i=l 

If  we  denote  the  left  side  of  (3.2.3)  by  G  (considered  as  a 
function  of  e  and  3)  and  write  GQ  as  the  partial  derivative 
of  G  with  respect  to  0  ,  then  if  3  is  Known  and  0  is  an 
estimate  for  0  ,  the  next  Newton-Raphson  recursive  estimate 
for  0  i s 


®2  V  G/Ge 


where  G  and  G  are  evaluated  at  0 

0  1 


In  the  same  way,  denoting  the  left  hand  side  of 

(3.2.4)  by  H,  if  0  and  3^  are  estimates  for  0  and  3  then 
the  next  estimates  are 


62  =  0J  +  (HGg  -  GHp)/(G0Hg  -  GgH0) 


6  =  6  +  ( GH  -  HG  ) / ( G  H  -  G  H J 

^2  0  0  03  30 


and 


. 
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where  again  the  right  hand  sides  are  evaluated  at  0 ^  . 

For  the  logistic  distribution  function  F  we  have 

H"  =  “ 3F ( x ) ( 1  - F ( x ) )  and  -  ( x- e ) F ( x ) ( 1  - F ( x  )  )  . 

d  p 

Using  these  we  can  obtain  the  quantities 

G  =  -  Y  w.F(X. ) ( 1-F(X.  )  ) 

0  L  i  i  i 

G3  =  l  wi ( X i - 0 ) F ( X i ) ( 1  - F ( X i  ) ) 

H  =  -  V  w.  X  .  F ( X . ) ( 1 -F ( X .  ) ) 

0  L  1  1  1  1 

fi  =  y  w. x .  (x.-e)F(x. )  ( i -f (x. ) ) 

3  L  1  1  T  1  1 

where  all  four  sums  are  taken  for  1  <  i  <  n.  The  same 

procedure  can  be  followed  with  unequal  starting  points  and 
increment  sizes.  The  equations  are  similar  with  the  w' s 
replaced  by  a  second  summation  sign. 

This  method  is  virtually  identical  to  Finney's  method. 
The  only  difference  is  that  the  sample  size  at  each  X^  is 
variable  and  the  samples  tend  to  be  closer  together  than  in 
the  classical  case.  Otherwise  the  maximum  likelihood 

equations  are  the  same  and  so  must  also  be  the  solutions. 
The  difference  in  the  iteration  scheme  is  due  to  the  fact 
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that  Finney  uses  a  first  term  expansion  of  the  distribution 
function  before  deriving  the  recursion. 

The  equations  of  course  are  much  simpler  for  the  one 
parameter  recursion  that  is  used  when  3  is  Known. 


1 


59 


3.3  Exponential  Estimates.  First  we  will  consider  the  case 
when  3  is  Known.  Theorem  2.2.1  shows  that  under  the  proper 
condi t ions 


P(m(e^D-1)  <  z)  ->  1  -  e"z 


so  for  this  section  we  will  assume  that 


(3.3.1)  P (eSD- 1  <  z)  =  1  -  e'mz. 

If  and  a  are  subscripted  by  i,  we  wi  1 1  write 


m  .  -&(Xi  -e) /,  &A . 

m.j  =  e  i.j  /  (e  i -  1  ) 


The  density  derived  from  (3.3.1)  is 


(3.3.2)  f(z)  =me 

For  K  such  variables  with  all  starting  points  equal  and  all 
increment  sizes  equal  (and  hence  only  one  common  value  for 
m)  we  have  the  joint  density 


(3.3.3)  L  =  f  (z  f  z 2,  .  .  zk) 


=  mk  exp-m  l  z  . . 

i  =  1  1 
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We  are  interested  in  estimates  of  ln(m).  Since  the  sum  of  K 
independent  Exponential  variables  with  mean  1/m  is  a  Gamma 
random  variable  with  mean  K/m  and  variance  K/m2 ,  we  may  use 
the  fact  that  the  expected  value  of  the  log  of  this  random 
variable  is 


where  ip  (K)  is  the  Digamma  function  which  for  an  integer  K 
is  equal  to 


y if  k= 1 


k-  1 

-  y +  l  1/i  if  k>1 . 


[Note:  y  =  Euler's  constant  =  .5772156645,..] 


Reca 1 1 i ng  that 


m  =  e-e(X1-6)/(eei-l) 


we  have  that 


E  [  1  n  X  Z-1  =  3  ( X  -e)  +  ln(e6A-1)  +  <Mk) 
i  =  1  1 


o  n 

and  so  if  we  put  D.  =  X  -  X  and  ep  i  - 1  =  z.  we  have  that 

i  l  n  • 


i 


' 
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(3.3.4)  6  =  X  ,  +  [ln(e  A  -  1  )  +  ^(K)  -  In  He  °i  -1)1/3 

1  i  =  1 


is  an  unbiased  estimate  of  Q  .  Since  3  is  a  function  of  the 
complete  sufficient  statistic,  it  follows  that  ^  is  a  UMVUE 
of  e . 

We  can  also  find  the  maximum  likelihood  estimate  for 
e.  Taking  the  log  of  (3.3.3)  and  using  subscripts  on  m  to 
allow  for  different  starting  points  and  increment  sizes 
gives 

k  k 

(3.3.5)  ln(L)  =  £  In  (mi  )  -  £  m-jZi. 

i  - 1  i  - 1 

Taking  the  partial  derivative  of  ln(L)  w.r.t.  0  gives 


1MikI=  mizi ' 


Setting  this  equal  to  0  and  solving  for  0  yields 


(3.3.6)  e  =  [ln(k)-ln  I  {e-B*!.  ( e  6°i  - 1 )  /  ( e  ^  1  - 1 )  >  ]  /  B 


This  estimate  of  course  can  be  used  even  with  unequal 
starting  points  and  increment  sizes.  However,  in  order  to 
compare  it  with  the  UMVUE,  let  us  assume  that  starting 
points  and  increment  sizes  are  the  same.  Then  (3.3.6) 


reduces  to 


■= 

l 
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(3.3.7)  0  =  [  ln(K)  -  In  £  {e  (e  3Di-  1 )  /  (e  3  A.  )  }  ]  / 

1  =  1 

=  X  +  [ln(eAA-1)  +  ln(k)-ln  J  (e^^i-1  ) 

1  i  =  l 

The  difference  between  the  UMVUE  and  the  MLE  is 

UMVUE  -  MLE  =  [y  (k)  -  1 n ( k ) ]/3 

which  is  a  contant  so  Var(UMVUE)  =  Var(MLE)  and  the  bias  of 
the  MLE  can  actually  be  calculated. 

When  3  is  not  known,  the  distribution  (3.3.1)  should 
be  rewritten 

(3.3.8)  P(D  <  d)  =  1  -  exp-m(e  - 1 ) 

and  the  corresponding  density  is 

(3.3.9)  f(d)  = 3  meAdexp-m(e^d -  1 ) 

so  the  new  log  likelihood  for  k  such  variables  is 

k  k  k 

(3.3.10)  ln(L)=kln(3)+  )’  ln(m.)  +3  J  d.  -  £  mi(e^di-1). 

i  - 1  1  i  =  l  1  i  =  1 

Since  we  already  have  the  MLE  for  0  (3.3.6)  we  need  only 

A 

solve  for  3.  (3.3.6)  gives  an  explicit  solution  for  0.  Thus 


' 
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when  we  take  the  partial  derivative  of  (3.3.10)  and  set  it 
equal  to  0,  we  may  substitute  (3.3.6)  into  the  equation 
leaving  it  a  function  only  of  3,  the  d' s,  A's,  and  X^'s. 
Thus  the  most  difficult  calculations  require  only  a  single 
parameter  recursion  to  solve  for  3.  This  value  may  then  be 
substituted  into  (3.3.6)  to  find  0.  Rather  than  give  these 
equations,  we  will  instead  consider  the  neater  solutions 
when  all  the  X  '  s  and  a's  are  equal.  In  this  case  (3.3.10) 


becomes 


k  k 

(3.3.11)  1  n  (  L  )  =  K 1  n  ( 3  )  +  Kln(m)  +  g  £  d  .  -  m  J  (e6d;.-1) 

i  =  1  i  =  l  1 


and  taking  the  partial  derivative  w.r.t.  3  yields 


(3.3.12)  -sjr1- 


=  k/s  +'  (k/m)|? 


i  =  1  i  =  1 


Substituting  (3.3.7)  into  m  gives 


=  k/  l  (eBdi-1) 


k 


and  substituting  this  into  (3.3.12)  gives 
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(3.3. 13) 


91  n  (  L  ) 

33 


k 

=  K/3  +  Id 
i  =  1  1 


k 


-  K(  l  die 
i  =  1 


3di 


)/  l  (e3di  -1). 

i  =  1 


Setting  this  equal  to  0  gives  an  equation  which  may  be 

/\ 

solved  for  3  numerically. 
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3.4  Extreme  Value  Estimates.  From  Corollary  2.4.1  we  Know 
that  under  certain  conditions 

P  (  3( X  N  -  0  )  +  In (e^A- 1 )  <z)  +  exp-e"  z . 

For  this  section  we  will  assume  that 


(3.4.1) 

P  (  X  ^  <  x)  =:  exp-exp-  3  (  x  -  {©-  [  ln(ef3A-1  )  ]/3>)  . 

Putting 

s  =  e  -  [ In ( e^A- 1 ) ] / 3 

we  can  reparametrize  (3.4.1)  to  the  form 


(3.4.2) 

P(>^  <  x)  =  exp-'exp- 3<x-s)  . 

This  is 

simply  an  Extreme  Value  distribution  (sometimes 

called 

a  Type  I  Extreme  Value  distribution)  with  mean  s  + 

if  /  3  and 

variance  jr2  /  (  6  3  2  )  .  The  estimates  we  are  interested 

in  are  the  moment  estimates,  found  by  solving 


(3.4.3) 

3  =  /  ( 6  times  the  sample  variance)2 

and 


(3.4.4) 

s  =  sample  mean  -  y/3- 

This  last  equation  can  be  solved  for  e  giving 


. 
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(3.4.5)  0  =  sample  mean  -  [y  -  ln(e  -1)]/3. 

We  will  not  consider  any  Maximum  Likelihood  estimates 
from  the  Extreme  Value  distribution  for  two  major  reasons. 
First,  the'  moment  estimates  for  1/3  and  s  given  above  are 
approximately  55%  and  95%  efficient  respectively  [12], 
While  this  may  be  low  as  far  as  estimating  3  is  concerned, 
the  high  efficiency  of  the  location  parameter  makes  it  very 
competetive,  since  the  bias  due  to  lack  of  fit  will 
dominate  the  effect  of  a  5%  loss  of  efficiency.  The  second 
reason  is  that  the  ML  estimates  from  the  Extreme  Value 
approximat ion  are  roughly  as  computationally  difficult  as 
those  from  the  Exponential  approximation,  and  since  the 
Exponential  approximation  is  better,  the  ML  estimates  based 
thereon  should  also  be  more  desirable.  It  is  primarily  the 
simplicity  of  the  moment  estimates  that  makes  them  so 
attractive.  We  will  see  that  they  also  perform  surprisingly 
well  under  certain  sampling  situations. 

Vie  could  of  course  calculate  estimates  from  the 
Truncated  Extreme  Value  approximation.  However,  since  this 
was  simply  a  transformation  of  the  Exponential 
approximation,  it  is  therefore  no  surprise  to  find  that  the 
UMVUE  of  e  when  3  is  known  and  the  ML  estimates  in  the 
general  case  are  identical  to  those  from  the  Exponential 
approximation,  and  we  will  not  reproduce  them  here. 
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CHAPTER  IV 


Relative  Performance  of  Estimates 

4.1  Introduct ion .  In  Chapter  III  we  presented  a  rather 
large  list  of  possible  estimates  which  could  be  used.  We 
have  actually  looked  at  several  others  but  because  of  their 
poor  overall  performance  we  did  not  feel  they  were  worth 
including  in  the  development  or  comparisons.  The  estimates 
can  be  classified  into  two  natural  groupings;  estimates  of 
6  or  r  when  3  is  Known  and  estimates  of  0 ,  rf  and  3  in  the 
general  case.  Section  4.2  will  treat  the  first  group  and 
section  4.3  will  concern  itself  with  the  second  group.  In 
addition  there  is  the  natural  question  of  what  would  happen 
if  the  model  were  incorrect.  That  is,  what  if  instead  of 
having  the  assumed  logit  model,  the  actual  underlying 
response  function  were  probit?  Such  a  question  is  obviously 
important  and  will  be  considered  in  section  4.4.  We  will 
then  end  the  chapter  with  a  comparison  of  First  Zero 
estimates  and  classical  logit  analysis  estimates  in  one 
particular  testing  situation. 

The  nature  of  the  estimates  generally  made  it 
impossible  to  get  the  explicit  moments  of  the  sample 
estimates.  As  a  result  most  of  the  comparisons  that  were 
made  were  done  by  means  of  Monte-Carlo,  simulations.  A 
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pseudo- random  number  generator,  the  one  supplied  by  the  APL 
installation  at  the  University  of  Alberta,  was  used  to 
simulate  the  large  samples  needed  to  be  able  to  compare  the 
bias  and  variance  of  the  different  estimates.  Wherever 
possible  we  used  paired  comparisons  and  we  also  used 
different  starting  seeds  on  different  runs  in  order  to  make 
the  results  of  separate  simulations  "independent".  We  also 
kept  track  of  the  starting  seeds  for  each  simulation  with 
the  result  that  each  of  the  simulations  is  completely 
reproducible. 

For  the  sake  of  standardizing  the  comparison 
procedure,  we  will  always  have  the  starting  points  equal 
and  increment  sizes  equal.  We  do  this  for  two  reasons. 
First,  having  a  randomly  varying  starting  point  and 
increment  size  would  add  variance  to  the  estimates  that 
might  obscure  relative  performance.  Second,  determining  an 
optimal  method  of  varying  X ^  and  A  sequentially  is  not  a 
simple  process,  and  very  likely  the  method  itself  would 
depend  on  the  type  of  estimate.  This  consideration  will  be 
discussed  again  in  the  next  chapter  in  the  sections  dealing 
with  recommendations  and  further  research. 

The  underlying  response  function  will  be  the  standard 
logistic  function 


F  (  x)  =  1/(1  +  e‘x  ) 


■ 
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which  has  e  and  3  equal  to  0  and  1  respectively.  Vie  tested 
for  three  different  starting  points:  LD75,  LD95,  and  LD99 
whose  numeric  values  are  In  3,  In  19,  and  In  99.  Three 
increment  sizes  were  used  as  well:  .02,  .1,  and  .5.  With 
shorter  programs  we  also  used  three  different  sample  sizes 
5,  10,  and  20.  For  longer  programs  fewer  sample  sizes  were 
used.  Here  sample  size  refers  to  the  number  of  O' s  in  the 
sample  and  not  the  actual  number  of  observations. 


, 
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4.2  Estimates  when  3  is  Known.  When  3  is  Known,  it  is 
irrelevant  which  root  is  to  be  estimated.  We  also  have  the 
atypical  statistical  situation  of  having  no  penalty  for 
extrapolation.  This  follows  since  if  r  is  the  p-fractile 
root  ( ie.  F ( r ) =p) ,  then 

r  =  e  +  { 1 n ( p/ [ 1 -p  3 ) } / 3 . 

Hence  if  e  is  any  estimate  of  e  then 

r  =  e  +  { ln(p/ [ 1 -p] ) }/ 3 

is  an  estimate  of  r  with 

Bias(r)  =  Biaslg)  and  Var(r)  =  Va n ( q ) . 

Thus  we  may  talk  simply  of  the  bias  of  an  estimate  without 
referring  to  the  root  it  is  estimating. 

Table  4.2.1  summarizes  the  results  of  twenty-seven 
separate  simulations  -  one  simulation  for  each  of  the 
starting  points,  increment  sizes,  and  sample  sizes.  Four 
estimates  were  considered:  the  Exact  MLE,  the  Exponential 
MLE,  the  Exponential  U.MVUE ,  and  the  Extreme  Value  Moment 
Estimate.  For  each  simulation,  one  thousand  samples  were 
generated  and  the  four  estimates  calculated  for  that 
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TABLE  4.2.1  COMPARISON  OF  BIAS  AND  VARIANCE  (IN  BRACKETS)  OF  ESTIMATES  WHEN  BETA  IS  KNOWN. 
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sample.  From  the  table  it  is  apparent  that  the  Exact  and 
Exponential  MLE' s  are  essentially  the  same  with  regard  to 
both  absolute  bias  and  variance,  except  perhaps  when  A  -js 
large,  in  which  case  the  Exact  MLE  might  have  slightly 
smaller  variance.  The  Exponential  UMVUE  has  the  same 
variance  as  the  Exponential  MLE  and  it  has  less  bias  than 
the  other  two.  The  Extreme  Value  Moment  Estimate  is  never 
really  competitive,  since  the  only  instances  that  its  bias 
is  acceptably  small  its  variance  is  as  high  as  that  of  the 
other  estimates.  The  Extreme  Value  and  Exponential 
estimates  are  all  very  straighforward  to  calculate 
explicitly,  while  the  Exact  estimates  require  a  recursion 
to  locate  the  root  of  the  ML  equation.  As  a  result,  the 
Exponential  UMVUE  seems  to  have  the  edge  as  far  as 
desirability  is  concerned,  but  it  may  only  be  used  when  all 
starting  points  and  increment  sizes  are  equal.  Otherwise, 
for  a  general  testing  situation  where  for  one  reason  or 
another  the  starting  points  or  increment  sizes  are  changed, 
the  Exponential  MLE  would  seem  to  be  the  best  choice.  The 
one  exception  to  this  might  be  in  the  situation  where  large 
values  of  A  were  to  be  used,  and  then  the  slight  decrease 
in  variance  might  justify  the  extra  computation  necessary 
to  use  the  Exact  MLE. 
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4.3  Estimates  when  r  is  not  Known.  In  the  general  case,  we 
will  only  consider  three  estimates:  the  Exact  (VILE ,  the 
Exponential  MLE,  and  the  Extreme  Value  Moment  Estimates. 
The  Exponential  UMVUE  may  only  be  used  when  3  is  Known. 
Four  quantities  must  be  taken  into  account  when  comparing 
different  estimates:  the  starting  point,  increment  size, 
sample  size,  and  also  the  root  to  be  estimated  since 
different  estimates  will  have  different  abilities  to 
extrapolate.  A  fifth  factor  which  may  be  quite  critical  in 
practice  is  the  computational  aspects  of  the  estimate.  The 
Exact  estimates  require  a  two  parameter  recursion  to  solve 
for  the  estimates,  while  the  Exponential  estimates  require 
only  a  one  parameter  recursion  and,  simplest  of  all,  the 
Extreme  Value  estimates  may  be  calculated  explicitly.  For 
this  reason,  the  Exact  estimates  are  the  most  sensitive  to 
choice  of  starting  approximations.  As  a  case  in  point,  for 
one  simulation  with  starting  point  LD95,  increment  size 
.02,  and  sample  size  10  we  used  the  starting  approximations 
of  0  and  1  (the  true  values)  for  0  and  3  respectively.  Even 
so,  we  had  to  generate  2369  samples  before  we  found  1000 
for  which  the  procedure  converged!  The  Exponential 
estimates  are  much  more  likely  to  yield  a  result.  In  the 
above  example,  1000  of  the  first  1115  samples  yielded 
solutions  for  the  Exponential  ML  estimates.  What  we  finally 
ended  up  doing  was  to  calculate  the  Exponential  estimates 
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for  the  sample  and  then  use  these  values  as  the  starting 
approximations  for  the  Exact  estimate  procedure.  When  this 
was  done,  the  proportion  of  samples  that  yielded  Exact 
estimates  increased  dramatically.  We  did  not  find  a  single 
case  when  the  Exact  procedure  converged  and  the  Exponential 
procedure  diverged.  The  problem  of  divergence  generally 
decreased  as  A  and  increased.  The  Extreme  Value 

estimates  were  undefined  only  when  all  the  stopping  points 
were  equal,  a  rare  occurrence  especially  when  the  sample 
sizes  were  at  least  10.  We  arbitrarily  decided  that  a 
procedure  would  be  considered  to  have  diverged  when  the 

/s 

value  of  3  had  gone  below  .1  or  above  10.  We  only  compared 
the  estimates  when  all  three  converged.  This  tends  to  favor 
firstly  the  Exact  estimates  and  secondly  the  Exponential 
estimates  since  they  were  not  penalized  in  any  way  for 

diverging.  We  had  considered  giving  the  limit  values  of 

/\ 

either  .1  or  10  to  3  when  the  procedure  diverged,  and  this 
would  have  had  the  effect  of  reversing  the  advantages  with 
the  Extreme  Value  coming  first  and  the  Exponential  coming 
second.  This  question  of  divergence  will  be  mentioned  again 
in  the  next  chapter  when  the  question  of  how  this  affects 
First  Zero  techniques  with  respect  to  Classical  methods 
under  cost  restrictions  is  discussed. 

Another  point  we  should  mention  is  that  of  adding  a 
Continuity  Correction  Factor  (CCF)  to  the  values  of  D  in 
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the  Exponential  and  Extreme  Value  approximations.  We  found 
that  a  CCF  did  not  help  either  the  fit  or  the  estimates  in 
the  Extreme  Value  case  and  so  we  did  not  include  one  there. 
Using  a  CCF  in  the  Exponential  estimate  for  r  when  3  was 
Known  actually  hurt  their  performance  and  again  we  did  not 
include  one.  However,  in  the  general  case  it  becomes 
apparent  that  leaving  out  the  CCF  adversely  affects  the 
performance  of  the  estimates.  This  is  primarily  due  to  the 
procedure  for  estimating  3  ,  which  is  why  the  CCF  was  not 
necessary  for  the  case  when  3  was  known.  Adding  the  CCF 
slightly  increases  the  bias  of  the  estimate  of  3  and 
decreases  its  variance.  These  two  factors  combine  to 
significantly  reduce  the  variance  of  the  estimates  of  the 
root  r.  For  these  reasons,  the  Exponential  MLE  will  be 
calculated  using  the  corrected  values.  To  each  value  of  D 
we  add  the  quantity  of  a/2  before  the  Exponential  estimates 
are  calculated. 

Tables  4.3.1  and  4.3.2  summarize  the  results  of 
eighteen  simulated  sampling  situations.  In  each  simulation 
samples  were  generated  until  1000  were  found  for  which  all 
three  estimates  converged.  For  each  sample,  all  three 
estimates  were  calculated  and  also  the  difference  between 
the  Exact  and  Exponential  MLE's.  Table  4.3.1  contains  the 
results  for  samples  containing  10  zeros  and  4.3.2  for 
samples  containing  20  zeros.  For  each  sample  estimates  were 
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TABLE  4.3.1  COMPARISON  OF  BIAS  AND  VARIANCE  (IN  BRACKETS)  OF  ESTIMATES  WHEN  BETA  IS  NOT  KNOWN 
AND  THE  NUMBER  OF  ZEROS  IN  EACH  SAMPLE  IS  10 
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TABLE  4.3.2  COMPARISON  OF  BIAS  AND  VARIANCE  (IN  BRACKETS)  OF  ESTIMATES  WHEN  BETA  IS  NOT  KNOWN 
AND  THE  NUMBER  OF  ZEROS  IN  EACH  SAMPLE  IS  20 
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calculated  for  3  and  the  two  closest  roots  to  the  sample 
points,  these  being  the  starting  point  and  the  next 
smallest  root.  In  general  we  found  that  extrapolation  very 
far  beyond  the  sampling  region  yielded  high  variances  and 
they  have  not  been  included  in  the  tables.  When  the 
starting  point  is  either  LD75  or  LD95  and  A  =  .02,  the 
sampling  region  stays  very  close  to  the  starting  point  and 
in  these  cases  we  only  gave  estimates  for  these  roots. 

From  the  tables,  the  following  general  remarks  may  be 

made : 

(i)  In  all  cases,  the  Exponential  estimate  of  3  is  better 
than  the  Exact.  It  has  both  less  absolute  bias  and  smaller 
var i ance . 

(ii)  The  Extreme  Value  estimate  of  the  root  has  less  Mean 
Square  Error  (MSE)  than  either  of  the  other  estimates  in 
all  estimates  based  on  samples  containing  10  zeros,  and 
only  loses  out  for  20  zeros  when  X  ^  =  LD99  and  A  =  .02.  It 
also  has  less  bias  when  X^  =  LD99  and  a  =  .5  for  20  zeros. 

(iii)  For  a  =  .02  and  X^  >  LD95  the  Exponential  and  Exact 
estimates  are  virtually  indistinguishable  with  the 
Exponential  perhaps  having  a  slight  edge.  The  results  are 
similar  for  X^  =  LD99  and  A  =  . 1. 

(iv)  As  the  number  of  zeros  in  the  sample  increases,  the 
bias  of  the  Exact  MLE  must  of  course  converge  to  zero.  It 
is  apparent  that  the  Exponential  MLE  will  also  converge  to 
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a  value  very  close  to  zero  in  most  cases  since  the 
approximation  is  very  good  for  X  >  LD95  and  A  <  . 1. 
However  the  bias  of  the  Extreme  Value  does  not  decrease  as 
the  number  of  zeros  increases,  but  instead  remains 
essentially  constant.  Thus  we  could  actually  calculate 
under  a  fixed  sampling  situation  for  what  number  of  zeros 
the  MSE  of  the  Exact  or  Exponential  estimates  would  become 
less  than  that  of  the  Extreme  Value.  However,  for  X^  large 
(say  at  least  LD95)  and  a  large  (eg.  .5)  this  would  require 
a  very  large  sample  before  the  ML  estimates  achieved  this 
goal  . 

Although  it  is  not  given  in  the  tables,  solutions  for 
both  the  Exact  and  Exponential  estimates  are  not  likely  to 
be  found  when  Xj  =  LD75  and  A  =  .02.  In  our  simulation  with 
10  zeros  in  the  sample,  2366  samples  were  generated  before 
1000  were  found  that  yielded  solutions  for  both  estimation 
procedures.  When  the  number  of  zeros  in  the  sample  was 
increased  to  20,  1000  of  the  first  1649  samples  gave 
convergent  solutions.  The  problem  diminished  as  X^  and/or  A 
and/or  the  number  of  zeros  increased.  The  reason  for  this 
can  be  seen  since  when  X^  is  central  and  A  is  small,  the 
samples  tend  to  get  selected  from  a  very  small  region.  This 
is  asymptotically  optimum  for  estimating  the  root  but  for 
smaller  sample  sizes  it  allows  considerable  variation  in 
the  estimate  of  3.  Recall  that  we  set  arbitrary  limits  on 
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the  permissible  values  of  3  at  .1  and  10.  For  such  a  small 
sampling  range  the  probability  of  3  falling  outside  this 
range  becomes  significantly  large.  Increasing  and/or  £ 
increases  the  sampling  range  and  reduces  the  problem. 
Increasing  the  sample  size  reduces  the  variance  of  the 
estimates  and  therefore  reduces  the  probability  of 
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4.4  Performance  of  Estimates  in  Probit  Analysis.  When  the 
response  function  is  not  logistic,  it  is  important  to  Know 
how  the  estimates  of  the  parameters  will  be  affected  by 
applying  logit  analysis  to  a  probit  model.  For  our 
simulations  using  the  logistic  distribution  we  used  the 
model  with  0  equal  to  0  and  3  equal  to  1.  This  produced  a 
distribution  with  mean  0  and  variance  tf 2 / 6  (=  1.6449...). 
For  our  Normal  response  function  we  will  use  a  Normal 
distribution  with  this  same  mean  and  variance  so  that  the 
biases,  variances,  and  mean  squared  error  of  the  estimates 
will  be  in  a  sense  comparable  to  those  in  the  tables  of  the 
previous  section.  The  increment  sizes  will  remain  the  same 
but  of  course  the  75th,  95th,  and  99th  percentiles  change. 
Their  values  for  the  above  Normal  distribution  are  ,865, 
2.11,  and  2.98  respectively.  These  values  are  much  closer 
to  zero  than  were  those  from  the  logistic  distribution,  due 
of  course  to  the  lighter  tails  of  the  Normal  distribution. 
An  alternate  method  of  comparing  the  two  response  functions 
would  be  to  choose  a  Normal  distribution  with  mean  0  and 
the  same  95th  percentile  as  the  logistic's.  If  the  prime 
root  of  interest  were  Known,  this  would  perhaps  be  a  better 
method,  but  since  we  were  interested  in  a  range  of  values, 
we  decided  to  set  the  means  and  variances  equal. 

As  in  the  last  section,  we  estimated  two  of  LD50, 
LD75,  LD95,  and  LD99,  depending  on  the  starting  point.  We 


< 

ft*] 

■ 


82 


have  not  performed  a  full  analysis,  but  have  only  included 
those  combinations  of  starting  points  and  increment  sizes 
for  which  one  of  the  Exponential  or  Extreme  Value  estimates 
performed  well  enough  to  be  used  in  actual  practice. 
Because  of  the  similarity  of  the  Exponential  and  Exact 
estimates,  as  wel 1  as  the  computational  difficulty  of  the 
latter,  we  have  omitted  it  in  the  comparison. 

Table  4.4.1  summarizes  the  results  of  the  simulation. 
From  the  table,  and  comparing  correspondi ng  values  from 
Table  4.3.1,  several  observations  may  be  made: 

(i)  The  bias  of  the  estimates  is  surprisingly  small 
considering  the  difference  in  the  models.  Of  course  these 
biases  wi 1 1  not  approach  0  as  the  sample  size  increases  as 
they  would  in  the  logit  model,  but  in  many  cases  the  number 
of  zeros  in  the  sample  must  be  very  large  before  the  square 
of  the  bias  becomes  a  major  factor  in  the  Mean  Squared 
Error  ( MSE ) . 

(ii)  The  variance  of  the  probit  estimates  is  generally 
smaller  than  for  the  logit  estimates.  This  can  be 
attributed  to  the  shape  of  the  curves  and  the  fact  that  the 
extreme  quantiles  of  this  Normal  distribution  are  much 
closer  together  than  the  logistic  extremes. 

(iii)  The  MSE  of  the  Extreme  Value  estimates  is  less  than 
the  MSE  of  the  Exponential  estimates  in  all  cases.  This  is 
due  to  the  lower  variance  of  the  EV  estimates.  However,  in 
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TABLE  4.4.1  BIAS  AND  VARIANCE  (IN  BRACKETS)  OF  EXPONENTIAL  AND 
EXTREME  VALUE  ESTIMATES  WHEN  THE  ACTUAL  UNDERLYING 
DISTRIBUTION  IS  NORMAL.  SAMPLES  CONTAIN  10  ZEROS. 
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this  case  the  Exponential  estimates  do  not  have  the 
advantage  of  smaller  asymptotic  bias  and  in  certain 
sampling  situations  it  is  likely  that  the  MSE  of  the 
Extreme  Value  estimates  would  always  be  less  than  that  of 
the  Exponential  estimates. 

(iv)  The  bias  is  a  function  of  the  mean  and  variance  of  the 
response  function  and  also  X  ,  a  *  and  the  root  to  be 
estimated.  Thus,  by  changing  X  ^  and  a appropr i a te ly  the 
estimates  could  be  made  to  be  unbiased  in  the  probit  model 
without  losing  their  asymptotic  properties  in  the  logit 
model.  We  will  mention  this  fact  again  in  the  next  chapter 
in  the  section  discussing  further  possible  research. 
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4.5  Comparison  of  First  Zero  and  Classical  Estimates.  It  is 


necessary  to  Know  not  only  how  the  First  Zero  estimates 
perform  in  relation  to  one  another,  but  also  how  they 
perform  with  respect  to  the  classical  analysis  estimates. 
This  could  be  thought  of  as  comparing  the  sampling 
procedures,  since  the  Exact  MLE  is  theoretically  the 
classical  estimate  adjusted  for  the  First  Zero  sampling 
technique.  To  do  this  we  ran  three  simulations  according  to 
reasonably  standard  procedures  as  recommended  by  Finney  in 
his  book  "Probit  Analysis".  We  sampled  at  four  different 
percentiles  in  each  simulation  and  took  equal  size  samples 
at  each  of  the  four  doses.  However ,  comparing  the  two  types 
of  estimates  does  pose  some  problems.  Classical  analysis 
has  fixed  sample  size  and  a  variable  number  of  zeros  in 
each  sample.  First  Zero  analysis  has  a  fixed  number  of 
zeros  and  a  variable  sample  size.  We  decided  to  compare 
samples  where  the  number  of  zeros  in  the  First  Zero  sample 
was  set  equal  to  the  expected  number  of  zeros  in  the 
Classical  sample,  modelling  the  situation  where  the  cost  of 
a  death  is  a  major  component  of  the  loss  function.  Also 
necessary  for  comparison  purposes  is  the  expected  First 
Zero  sample  size  and  this  is  also  provided  in  the  table. 
Both  sampling  methods  are  most  effective  when  the  sampling 
points  bracket  the  point  to  be  estimated.  For  this  reason 
we  restricted  our  attention  to  the  case  of  estimating  LD95 


■ 
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and  starting  at  LD99  for  the  First  Zero  estimates.  The 
First  Zero  values  are  reproduced  here  again  from  Table 
4.3.2  to  facilitate  the  comparison.  We  considered  only 
samples  containing  20  zeros  in  the  First  Zero  analysis  and 
samples  having  a  mean  of  20  zeros  in  the  classical 
analysis.  We  also  included  the  estimates  of  3  in  the 
comparison,  making  3  and  LD95  the  parameters  of  interest. 

The  results  are  summarized  in  Table  4.5.1.  From  the 
table,  several  comments  may  be  made.  Some  of  these  are: 

(i)  For  a  fixed  number  of  zeros  and  approximately  the  same 
sample  size,  each  of  the  First  Zero  estimates  of  both  3  and 
LD95  are  better  than  the  correspondi ng  Classical  estimates. 
The  only  exception  to  this  is  the  Extreme  Value  estimate  of 
3  when  A  =  . 1 . 

(ii)  Since  the  Exact  ML  estimates  are  essentially  the  same 
as  the  Classical  ones,  it  would  seem  that  that  the  First 
Zero  procedure  might  be  a  more  efficient  sampling  technique 
for  estimating  extreme  quantiles  in  logit  analysis,  at 
least  when  the  number  of  zeros  is  a  limiting  factor  in 
determining  the  sampling  method. 

(iii)  It  appears  that  under  certain  sampling  situations, 
the  Extreme  Value  moment  estimates  can  provide  a 
significant  increase  in  efficiency  over  the  usual  Maximum 
Likelihood  methods  for  even  reasonably  large  samples.  By 
the  same  token,  in  many  instances  the  Exponential  estimates 
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TABLE  4.5.1 
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BIAS  AND  VARIANCE  (IN  BRACKETS)  OF  ESTIMATES 
OF  B  AND  LD95 .  FIRST  ZERO  SAMPLES  CONTAIN  20 
ZEROS.  CLASSICAL  ANALYSIS  SAMPLES  HAVE  A  MEAN 
OF  20  ZEROS  PER  SAMPLE. 


FIRST  ZERO  ESTIMATES 


TYPE  OF 
ESTIMATE 

S 

A  ' 

ESTIMATE 

OF  3 

ESTIMATE 
OF  LD95 

MEAN  #  OBS 
PER  SAMPLE 

EXACT  MLE 

LD99 

.  1 

.  105 

( .088) 

.007 
( .060) 

424.6 

EXP  MLE 

LD99 

.  1 

.061 
(  .078) 

-  .031 
(  .059) 

424.6 

EV  MOM 

LD99 

.  1 

.405 
( .049) 

-.057 
(  .043) 

424.6 

EXACT  MLE 

LD99 

.5 

.092 
( .062) 

-  .050 
( . 186) 

161  .0 

EXP  MLE 

LD99 

.  5 

-  .059 
( .040) 

-  .026 
( . 188) 

161  .0 

EV  MOM 

LD99 

.5 

.072 
( .034) 

.012 

(.111) 

161  .0 

CLASSICAL 

ANALYSIS  ESTIMATES 

SAMPLING  METHOD 

ESTIMATE 

OF  3 

ESTIMATE 
OF  LD95 

#  OF  OBS 
PER  SAMPLE 

100  OBS  AT 

LD92 ,94,96 

,98 

.119 
(  .256) 

-  .053 
( . 107) 

400 

50  OBS  AT 

LD8 1 ,87,93, 

99 

.  147 
(  .223) 

-  .031 
( . 195) 

200 

10  OBS  AT 

LD20, 40,60, 

80 

.  127 
( . 194) 

.  140 
(2.  14) 

40 

' 
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are  as  good  or  better  than  Exact  or  Classical  estimates 
with  the  added  bonus  of  being  easier  to  calculate.  The 
Extreme  Value  estimates  are  by  far  the  most  attractive 
computationally,  however. 
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CHAPTER  V 


Conclusions 

5.1  Cost,  Loss,  and  Risk  Considerations.  The  First  Zero 

sampling  technique  is  desirable  for  two  reasons:  sampling 

is  done  near  the  root  of  interest  and  the  number  of  zeros 

in  the  sample  can  be  controlled.  As  mentioned  previously, 

by  changing  the  values  of  and  A,  the  proportion  of  zeros 

to  observations  in  the  sample  and  the  MSE  of  the  estimates 

can  be  varied.  This  becomes  important  when  the  relative 

costs  of  taking  an  observation  and  having  a  zero  are 

considered.  As  an  example,  let  us  say  that  the  cost  per 

observation  is  C  and  the  additional  cost  of  a  zero  is  C^ 

(as  in  the  case  when  a  zero  corresponds  to  a  death).  If  the 

total  cost  of  the  experiment  was  to  be  C  and  the  number  of 

observations  were  N  and  the  number  of  zeros  were  N^,  then 

we  would  want  to  minimize  the  MSE  (or  perhaps  some  other 

loss  function)  subject  to 

N  C  +  N  C  <  C. 

11  2  2 

In  practice,  it  is  rare  that  both  N  and  N  would  be  known 

12 

before  the  experiment  begins;  generally  one  is  fixed  and 
the  other  is  variable.  Hence,  another  possibility  is  to 
minimize  the  MSE  subject  to  the  alternate  constraint 

E(N1)C1  +  E(N2 )C2  <  C 

where  E ( X )  is  the  expected  value  of  the  random  variable  X. 
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In  First  Zero  sampling  it  is  then  necessary  to  Know  the 
expected  number  of  observations  per  zero.  For  the  values  we 
have  discussed  these  quantities  have  been  calculated  and 
are  presented  in  Table  5.1.1.  For  classical  sampling,  the 
total  number  of  observations  would  be  set  and  the  expected 
number  of  zeros  calculated. 

Let  us  illustrate  the  above  by  means  of  a  rather 
extensive  example  using  the  latter  of  the  two  above 
constraint  conditions.  We  assume  we  wish  to  estimate  LD95. 
Further  we  assume  that  C^=1,  C2  =  5,  and  C=300,  We  cannot  do 
a  complete  analysis  of  the  optimal  sampling  strategy  for 
each  estimate,  but  we  shall  try  to  give  a  represen tat i ve 
survey  of  the  competing  estimates.  To  calculate  the 
expected  cost,  we  had  to  fix  the  number  of  zeros  in  the 
First  Zero  samples  and  calculate  the  expected  number  of 
observations,  while  in  the  Classical  samples  the  number  of 
observations  was  fixed  and  the  expected  number  of  zeros 
calculated.  For  the  sampling  method  we  indicate  the 
starting  point  and  increment  size  for  First  Zero  sampling 
and  the  dose  levels  and  number  of  replicates  at  each  level 
for  the  Classical  technique.  The  results  are  summarized  in 
Table  5.1.2.  The  MSE  was  calculated  by  generating  1000 
samples  for  each  of  the  estimates.  The  simulations  used 
different  seeds  and  thus  may  be  considered  independent. 
Again  estimates  were  not  penalized  in  the  case  of 
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XPECTED 

NUMBER  OF 

OBSERVATIONS 

PER  ZERO  IN  THE 

LOGIT 

MODEL 

X1 

INCREMENT 

SIZE 

.02 

.  1 

.5 

LD75 

3.84 

3.42 

2.63 

LD95 

15.6 

10.0 

5.  12 

LD99 

47.0 

21.2 

8.05 

TABLE  5. 

1  .2 

COMPARISON  OF  MSE  OF  VARIOUS ‘ESTIMATES 
OF  LD95  WITH  COST  RESTRAINTS 


TYPE  OF 
ESTIMATE 

SAMPLING  METHOD 

E  ( N  ) 

e(n2) 

E ( COST  ) 

MSE 

EXPONENTIAL 

LD95 ,  .02 

218 

14 

288 

.259 

EXPONENTIAL 

LD99 ,  .1 

233 

1  1 

288 

.  126 

EXPONENTIAL 

LD99 ,  .5 

185 

23 

300 

.  159 

EXT  VALUE 

LD95 ,  .02 

218 

14 

288 

.076 

EXT  VALUE 

LD99 ,  .1 

233 

1  1 

288 

.092 

EXT  VALUE 

LD99 ,  .5 

185 

23 

300 

.099 

CLASSICAL 

60@LD92 ,94,96,98 

240 

12 

300 

.253 

CLASSICAL 

500LD81 ,87,93,99 

200 

20 

300 

.203 

CLASSICAL 

2 1 OLD20 ,40,60,80 

84 

42 

294 

.905 

CLASSICAL 

1 20@LD92 , 98 

240 

12 

300 

.  167 

CLASSICAL 

1 00@LD82 , 98 

200 

20 

300 

.  155 

CLASSICAL 

42GLD 10,90 

84 

42 

294 

.380 

i 

1 
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divergence,  since  divergent  solutions  were  ignored.  This 

question  is  important  enough  to  warrant  further  attention. 

Both  First  Zero  and  Classical  techniques  will  not  yield 

solutions  for  certain  samples.  The  only  general  statement 

that  can  be  made  is  that  the  Extreme  Value  estimates 

diverge  less  in  all  situations  in  Table  5.1.2  than  any  of 

the  other  estimates.  Between  the  Exponential  and  Classical 

estimates  it  depends  on  the  situation;  some  Exponential 

estimates  diverge  less  than  some  Classical  estimates  and 

# 

vice  versa.  The  bottom  line  perhaps  is  comparing  those 
cases  with  the  least  MSE  for  each  technique.  In  this  case 
the  Exponential  both  had  smaller  MSE  and  diverged  less  than 
the  Classical  estimates. 

From  the  table,  several  observations  may  be  made: 

(i)  The  Extreme  value,  estimate  is  clearly  the  winner  in 
this  situation.  It  is  also  remarkably  unaffected  by  choice 
of  starting  point  and  increment  size.  Although  the 
combination  giving  the  smallest  MSE  was  X^=LD95  and  a  =.02, 
the  major  part  of  the  MSE  was  the  bias  of  the  estimate. 
Thus  perhaps  one  of  the  other  combinations  would  be 
preferrable  in  a  larger  experiment. 

(ii)  Although  not  given  in  the  table,  it  should  be 
mentioned  that  the  bias  of  all  the  estimates  except  one  was 
quite  small.  The  absolute  bias  of  the  Exponential  estimates 
was  less  than  .025,  and  that  of  the  Classical  estimates 
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less  than  .1.  As  mentioned  above,  the  Extreme  Value 
estimate  for  =  L D 9 5  and  A =.02  was  the  worst  at  -.259,  but 
the  other  EV  estimates  had  absolute  biases  less  than  .05. 

(iii)  We  have  included  various  allocations  of  observations 
for  the  Classical  method  in  patterns  roughly  similar  to 
those  proposed  for  example  by  Finney,  Wetherill,  and 
Tsutakawa.  For  this  situation,  the  minimum  MSE  is  achieved 
by  placing  two  groups  of  100  observations  at  each  of  LD82 
and  LD98.  Apparently,  fewer  dose  levels  with  larger  numbers 
of  observations  per  dose  seems  to  be  more  efficient. 

(iv)  The  Exponential  estimate,  while  not  matching  the 
Extreme  Value  estimate,  seems  to  be  better  than  the 
Classical  estimates.  Since  both  are  Maximum  Likelihood 
estimates,  this  would  seem  to  indicate  that  First  Zero 
sampling  is  a  more  efficient  method  of  placing  observations 
than  the  Classical  one,  at  least  under  certain  restraint 
conditions.  An  obvious  modification  to  Classical  sampling 
would  be  to  vary  the  number  of  observations  at  each  dose. 
This  is  essentially  what  is  accomplished  by  Wetherill's 
UDTR  rule  and  a  comparison  of  the  performance  of  this 
estimate  to  First  Zero  estimates  for  extreme  quantiles 
would  be  an  interesting  topic  of  further  research. 

(v)  In  actual  practice,  the  true  value  of  the  sampling 
doses  and  spacing  will  not  be  known  but  will  be  randomly 
distributed.  Their  distribution  will  either  be  determined 


■ 
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by  prior  Knowledge,  or  by  previous  iterations  of  the  root 
finding  process.  These  random  variations  introduce  another 
component  into  the  MSE,  but  the  expected  MSE  is  just  the 
expected  value  of  the  MSE  given  these  values.  It  is 
therefore  important  that  the  MSE  not  change  much  if  the 
doses  or  increments  are  changed.  The  Extreme  Value  estimate 
has  this  property,  and  both  the  Exponential  and  Extreme 
Value  estimates  have  relatively  stable  cost  functions, 
whereas  the  Classical  method  does  not. 
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5.2  Sequential  Approximation  Procedures.  Our  original  aim 
when  we  started  research  in  this  area  was  to  develop  a 
sequential  approximation  procedure  that  improved  on  the 
Robbins  -  Monroe  one  defined  by  (1.3.1).  In  this  we  were 
not  entirely  successful  if  one  requires  that  the  procedure 
must  both  (a)  have  the  nth  dose  contain  all  the  information 
from  the  first  n-1  observations  and  (b)  use  the  nth  dose  as 
the  present  estimate  of  the  root  in  order  to  be  considered 
a  sequential  estimation  procedure.  Vie  do  not  believe  this 
goal  can  be  achieved  without  significant  loss  of  efficiency 
in  the  estimates.  The  reason  for  this  is  the  type  of  errors 
that  occur  at  extreme  quantiles.  For  example,  sampling  near 
LD95  will  result  in  two  types  of  errors,  the  one  resulting 
from  observing  a  one  which  gives  an  error  of  .05  with 
probability  .95  and  that  of  observing  a  zero  which  gives  an 
error  of  .95  with  probability  .05.  Thus  the  error 
associated  with  a  zero  is  nineteen  times  as  large  as  that 
associated  with  a  one.  Any  Robbins  -  Monro  type  procedure 
will  cause  a  jump  to  the  right  nineteen  times  as  large  as 
the  jump  to  the  left  when  a  zero  is  observed.  Hoping  to  be 
in  the  vicinity  of  the  root  under  these  circumstances  would 
require  a  very  small  increment  size. 

However,  our  methods  could  be  considered  "almost" 
sequential  approximation  procedures.  Anytime  after  the 
first  zero  is  observed  you  will  have  a  continuously  updated 
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estimate  of  the  root.  The  sampling  procedure  itself 
requires  sequential  sampling.  And  the  choice  of  will 
generally  be  made  in  such  a  way  that  the  sampling  will  be 
done  in  the  vicinity  of  the  root.  If  simplicity  is  a  prime 
consideration,  then  the  Extreme  Value  Moment  estimates 
should  be  given  consideration.  This  procedure  only  requires 
recording  the  sum  and  sum  of  squares  of  the  stopping 
points,  and  could  easily  be  programmed  into  a  pocket 
calculator  such  as  the  Texas  Instruments  model  57,  58,  or 
59.  If  rigor  is  required  and  computing  facilities  are 
available  then  the  Exact  ML  estimates  offer  full  efficiency 
with  perhaps  even  the  added  advantage  of  more  efficiency 
than  the  classical  analysis.  In  either  case  the  procedure 
could  be  programmed  to  indicate  what  dose  is  required  at 
each  stage  and  when  the  desired  accuracy  was  achieved.  The 
methods  are  also  fairly  resistant  to  incorrect  assumptions 
about  the  form  of  the  response  function  which  makes  them 
desirable  in  situations  where  there  are  doubts  about  the 


model . 


fi 


97 


5.3  Recommend at  ions .  Two  cases  should  be  considered  here:  3 
known  and  3  unknown.  When  3  is  known  the  best  choice  would 
be  the  Exponential  UMVUE  since  its  bias  is  generally  less 
than  the  other  two.  Starting  at  a  fairly  central  value  with 
a  smaller  increment  size  seems  to  give  smaller  variance 
than  starting  at  a  more  extreme  quantile  and  using  a  larger 
increment.  If  the  original  starting  point  is  i nappropr i ate 
X  |  may  be  changed  and  the  observations  used  in  the 
Exponential  ML  estimate  in  the  more  general  sampling 
situation. 

When  3  is  not  known,  the  choice  of  estimate  depends  on 
more  factors.  Already  mentioned  in  the  previous  section  is 
simplicity.  Here  the  Extreme  Value  estimate  wins  easily. 
The  choice  of  A  is  critical  since  it  will  govern  the  number 
of  observations  per  zero,  and  the  .  estimate  makes  no 
allowance  for  changing  A  .  The  original  starting  point  is 
not  a  major  concern,  however,  since  the  starting  point  is 
not  used  in  calculating  the  estimate.  X}  must  of  course 
eventually  be  set  to  a  sufficiently  large  value  to  ensure 
reasonably  small  bias. 

The  optimal  values  of  X^  and  A  will  depend  on  cost 
considerations  and  also  the  unknown  parameters  0  and  3.  In 
practice  these  parameters  will  not  be  known  and  so  X ^  and  A 
will  have  to  be  changed  during  the  sampling.  The  general 
technique  would  be  to  sample  until  a  0  is  observed.  Based 
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on  the  sampling  up  to  that  point  estimates  of  the  root  r 
and  3  would  be  calculated.  Using  these  estimates  the 
optimal  values  of  Xj  and  A  would  be  found  for  the  next  set 
of  observations.  Following  this  procedure  would  guarantee 
that  and  A  would  converge  to  their  optimal  values.  What 
the  optimal  values  would  be  is  an  open  question  which  we 
mention  in  the  next  section. 

Another  consideration  is  the  loss  function  to  be  used. 
If  the  number  of  zeros  is  small  to  medium  and  squared  error 
loss  is  used  then  again  the  Extreme  Value  Moment  estimates 
will  likely  be  the  best  choice.  However,  if  unbiasedness  is 
important,  then  one  of  the  other  two  estimates  will  be 
better,  especially  if  the  root  is  fairly  central.  If  the 
sample  is  small  then  a  large  starting  point  and  large 
increment  size  should  be  used  to  increase  the  probability 
of  obtaining  a  convergent  solution.  If  a  large  sample  is  to 
be  used,  then  a  smaller  increment  should  eventually  be  used 
in  order  to  concentrate  the  observations  about  the  root.  If 
3  is  to  be  estimated,  then  either  the  Exponential  MLE  (less 
bias)  or  the  Extreme  Value  Moment  estimate  (less  MSE  under 
certain  conditions)  should  be  used. 

Finally,  the  cost  considerations  must  be  taken  into 
account,  in  the  choice  of  Xj  and  A  as  outlined  in  Section 
5.1  and  these  factors  are  fixed  for  all  the  estimates.  What 
determines  the  choice  of  estimate  will  be  which  one 
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performs  best  (whatever  the  loss  function)  under  the  given 
constraints. 

' 
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5.4  Further  Research.  Here  we  would  like  to  mention  briefly 
some  questions  that  this  thesis  raises  and  possible 
research  problems  arising  from  them. 

(a)  In  the  case  when  3  is  known,  it  would  be  useful  to  know 
the  exact  relationship  between  the  MSE  and  the  choice  of 
and  A  .  This  would  require  straightforward  but  extensive  and 
large  scale  simulations  to  obtain  either  a  function  or 
graph  of  the  one  as  a  function  of  the  other  two. 

(b)  As  mentioned  previously,  it  would  be  useful  to  know  the 

optimal  values  of  and  A  for  estimating  the  root. 

Unfortunately  these  values  will  be  functions  of  the  root  to 
be  estimated,  the  loss  function  used,  the  distribution 
parameters,  and  the  cost  restraints.  It  would  be  impossible 
to  tabulate  such  values  under  all  possible  combinations  of 
these  factors  and  so  a  study  on  the  relationships  between 
these  factors  and  the  choice  of  X^  and  A  would  have  to  be 
made  in  order  to  understand  their  behavior  in  the  general 
case. 

(c)  Along  this  line  it  would  be  useful  to  distinguish 
between  models  while  sampling  only  from  the  extreme 
quantiles.  To  do  this  a  goodness  of  fit  test  could  be  used 
to  determine  whether  the  distribution  of  the  stopping  point 
(for  example)  was  approximately  Extreme  Value,  but  what  is 
the  limiting  First  Zero  distribution  in  the  probit  model? 
It  might  be  possible  to  do  such  a  test  using  many  fewer 
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zeros  than  are  now  required  if  the  limiting  distributions 
of  the  different  response  functions  were  Known.  The 
advantage  here  is  that  the  distributions  differ  the  most  in 
the  tails  and  that  is  precisely  where  the  First  Zero  method 
concentrates  the  observations. 

(d)  Finally  it  would  be  nice  to  Know  if  the  First  Zero 
sampling  method  is  more  efficient  than  classical 
techniques.  Under  the  requirements  that  the  number  of 
observat ions/zero  be  equal  to  the  expected  number  of 
observat ions/zero,  is  the  proposed  procedure  the  best  way 
of  selecting  a  sample  from  the  tails  of  a  response 
function?  Again,  this  would  require  a  large  scale 
simulation  involving  running  the  classical  analysis  with 
different  sample  sizes  at  each  dose  and  different  numbers 


of  doses . 
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APPENDIX 


Fop  completeness  we  include  here  the  program  listing  of  all 
the  APL  programs  used  in  the  calculations  and  simulations 
presented  in  the  thesis.  Below  is  an  index  to  the  programs. 
The  programs  are  listed  alphabetically,  one  per  page  unless 
more  than  one  page  is  required.  The  index  indicates  the 
programs  that  were  used  for  each  table  or  figure. 


Where  Used 


List  of  Programs  Used 


F igures  2.3. 1 -2 . 3 . 9 
F igures  2.5. 1 -2 . 5 . 9 
Table  4.2.1 
Table  4.3.1 
Table  4.3.2 
Table  4.4.1 
Table  4.5.1 
Table  5.1.1 
Table  5.1.2 


SINKEXP,  GEN 

SINKEV,  GEN 

THETA,  GEN 

BETA,  GEN,  BEXPML 

BETA,  GEN.  BEXPML 

NOREV,  NORSIM,  NORMAL,  GENOR 

FIN,  BETA,  GEN,  BEXPML 

GEN 

FIN,  EVMOM,  EVEX,  GEN,  BEXPML 
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[0] 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 

[29] 

[30] 

[31] 

[32] 

[33] 

[34] 

[35] 

[36] 

[37] 

[38] 

[39] 

[40] 

[41] 

[42] 

[43] 

[44] 

[45] 

[46] 


XI  BETA 
T<-UAIl  2] 

t  t 


DEL ;S;SS;T ; CNT ; XI ;  I ; NTH ; NB \T H ; B ; P ; PI ; FAC ; DEN ; 

/  ;  B  ;  C ;  D ;  E ;  NC ;  F55  \  R  \  F  \  G  \  H  \N  \WT  ;DD 


'THIS  FUNCTION  FINDS  THE  EXACT ,  EXPONENT  I A  L  ,  yJ/VD  EXTR 

EME ' 

'VALUE  ESTIMATES  FOR  BETA ,  THETA  y  LD  75,  L£>95y  LF9 

9  .  * 

'SPECIFICATIONS:  ' 

’  XI  DELTA  SAMPLE  SIZE  SIMULATION  SIZE 

SEED ' 

6372120150  20  0  vXl , DEL , NUMDEATHS , SIMSIZE 9QRL 
S<-SS+E<-  6  5  ptfC-'-CWT’^O 
P«-l  0  0  0  0  x  +  \GEN  XI  y  DEL 
R+®  1  3  19  99 

MAINLOOP: N+ 1  +  +  / ( ? NUMDE AT H S p  1  0  0  0  0  ) °  .  >P 
n  CALCULATE  THE  EXPONENTIAL  ESTIMATES 
B+-BEXPML  DD+-DE L*N -  0 . 5 

T H<-X  1  +  ( -s- B  )  x ®NU MDE A T HS x  (  "  1  +  *B*DEL  )  ±  +  /  ~  1  +  * B*DD 
NC^NC  +  ERR*-(B<0  .  1  )  v5>9  .  9 
-^ERR/ MAIN  LOOP 
EiUl  +  (-l  yR)+B  yTH  +  R+B 
p  CALCULATE  THE  EXACT  ESTIMATES 
I*-  0 

WT<-+/ (  X  r/AO  o  .  <// 

XJ^Xl+PFIxi- ipl/71 
LOOP:  +  {  10  >!<-!+ 1  )  / CONTINUE 
-+(1<NC+NC+1  )  /MAIN LOOP 
CONTINUE  :PI<-l  +  l  +  *-  B*XI-TH 
A+-+  /FAC+PIxF<-WTxl-PI 
H<-+  /  F  AC*  XI 
C<-  ’+/F  AC  xXI*XI 
D++/XIIN1 
G*-+/XI*F 
F<-NUMDE  AT  HS  -  +  /  F 

NC+NC+ERR+O  .  0001  >  \ DEN+ ( B*A*C - T H* H) - BxH* H - T H*A 
-+ERR/ MAIN  LOOP 

NT  H+-T  H+  {  +DEN )  x  (  ( H-TH*A  )  *G  -  D)  +F*C  -T  H*H 
NB+-B+  ( +DEN)  x(FxBxi/)+Bx^xC-£) 

NC+-NC  +  ERR+-  (  NB  <  0  .  1  )  v  ( NB  >  1  0  )  v  2  5  <  I  NTH 
-+ERR/MAINL00P 

->(  (0 . 01  <  I  ( B<-NB  )  -  £  )  v  0  .  0 1  <  I  (  TH+-NTH )  -  TH)  / LOOP 
F[3;  ]<-(-!  ,  R)  +B  ,  TH+R+B 
p  CALCULATE  THE  EXTREME  VALUE  ESTIMATES 
B<-(  ( NUMDE  AT  H  S*  +  i  DD*DD )  -  (  A  +  +  /  DD)  *2) +NUMDE  AT  H  S*  NUMDE  AT  II 

5-1 

£«-  (  Ol  )*(  6x5)  *0 . 5 

75«-(Xl  +  (  0  . 5* DEL)- A* NUMDE AT HS)+(±B) x“ 0 . 5772157+©- l  +  *£x 

DEL 

El5;']  +  (-lyR)+ByTH  +  R+B 
p  CALCULATE  THE  DIFFERENCES 
F[  2  ;  ]<-F[  1  ;  ]  -  F[  3  ;  ] 

F[4;>F[3;]-F[5;] 
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[47] 

[48] 

[49] 

[50] 

[51] 

[52] 

[53] 

[54] 

[55] 

[56] 


£[6;  ]«-£[  5;  1-Ell  ;  ] 

S+-S+E 

SS<~SS+E*E 

->  (  SIM  SI  ZE>CNT<~CNT+1  )  /MAINLOOP 
’  BIAS 

VARIANCE' 

100p’  BETA  THETA  LD  75  L£>95  LD99 

t 

A  +  6  8  p  *  EX PON  EXACT  EX  VAL 

» 

8  3  vS+N^SIMSIZE 
A  ,  (S  8  p»  »),  8  3  ▼(SS*W-1)- 

'CW  TIM#  ’  ,  (TO  .001xDyU[2]-n ,  (10  0  vNC),'  DIVERGED ' 


* 

\ 
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[0] 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 
[7] 
C8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 
[19] 


B+-BEXPML  D\K\  EBD ;  SEBD  ;DBAR;FB;  FPB  ;  NEWB  ;  DEBD 
r THIS  FUNCTION  USES  A  NEWT  ON  -  RAPE SO N  ITERATION 
n TO  FIND  THE  ML  SOLUTION  FOR  BETA  FROM  THE 
r EXPONENTIAL  APPROXIMATION .  SINCE  THE  VALUE 
*OF  BET A-0  IS  A  FALSE  SOLUTION ,  A  PRELIMINARY 
r SEARCH  IS  MADE  TO  FIND  A  REASONABLE  STARTING 
r POINT,  IF  NO  SUCH  POINT  IS  FOUND ,  THE  VALUE 
r OF  BETA< 0.1  OR  BETA =10  IS  RETURNED . 

K+-pD 
B+ 0 

AG  A  IN:B<-B  +  2 
-*(£>9 . 9  )  /0 

LOOP:  SEBD+-+  / EBD<-*B*D 
DEBD+-+  /  D*EBD 
DBAR+-  (  +  / D)  ±K 

FB*~  ( SEBD-K )  +  (BxDBARxSEBD-K )  -B*DEBD 
-*•  (  0 <FB  )  /  AGAIN 

FPB*-(DBAR*  SEBD-K)  +  (B*DBARxDEBD)  -  B*  +  /  EBD*D*2 
NEWB+B-FBVFPB 

->•((£>  0  .  1  )  a  o  .  0 1  <  I  (B+-NEWB)  -B)  /  LOOP 
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[0] 

[1] 

[2] 

[3] 

m 

[5] 

[6] 

[7] 

[8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 

[29] 

[30] 

[31] 

[32] 

[33] 

[34] 

[35] 

[36] 

[37] 

[38] 

[39] 

[40] 

[41] 

[42] 

[43] 


EVEX ;  7  ;  A'  1  ; DEL ;ND;NS;S;SS ; CUT ; P ; 1 ; T H ; B ; N ; B ; D ; E ; 1C ; ERR ; 


R\AV\VAR\SKIP\A 

T+047[2] 

'THIS  FUNCTION  COMPARES  THE  MAXIMUM  LIKELIHOOD  E ST  IMA 

TES' 

' FROM  THE  EXPONENTIAL  APPROXIMATION  TO  THE  MOMENT  EST 

IMATES' 

'FROM  THE  EXTREME  VALUE  APPROXIMATION . ' 

' ENTER  THE  STARTING  POINT  XI:' 

Xl-D 

'ENTER  THE  INCREMENT  SIZE  DELTA:' 


DEL--  □ 


'ENTER  THE  NUMBER  OF  DEATHS  IN  EACH  SAMPLE :' 

ND+-D 

' ENTER  THE  NUMBER  OF  SAMPLES  IN  THE  SIMULATION: ' 

NS+-U 

'ENTER  THE  RANDOM  NUMBER  GENERATING  SEED: ' 

BRL+-U 

' ENTER  A  1  IF  YOU  DO  NOT  WISH  TO  SEE  INTERMEDIATE  RES 

ULTS , ’ 

*  A  OIF  TOU  DO: ' 

SKIP+ □ 

S+SS+E+  3  5  p  ZC--Z  pCNT<-0 
P<-1  0000x  +  \GEN  XI, DEL 
R+®  1  3  19  99 

AGAIN:  D<-DEL*  0  .  5  +  +  /  (  ?/7Dpl  00  00  )  °  .  >P 
B+-BEXPML  D 

T H+-X 1  +  (  tB  )  x©Arz?x  (  "  1  +  *B*DEL  )*+/-  1  -  *£x£ 


ZCll  3  +  ZC[  1  ]  +ERR+-(B<C  .  1  )  vB>9 . 9 
E  [1  ;  ]<-(  1  -ERR )  * B  ,  T H  +  R+  B 

VAR<-0 . 01  [  (  (+/D*D)*ND-  1  )  -  (  (  A  (  +  ID  )  RND )  *  2  )  *NDRND-  1 

ZC[2]+-ZC[2  ]+ERR+-(B<C  .  1  )  v  9  .  9  <5- ( Ol  )  -f  (  6  x  VAR )  *  0 . 5 

El  2  ;  ]<-(  1  -ERR)  *B  ,  (XI  -  AV-DEL  +  2  )  +  (-fB)xA’+“0.577216  +  ®"l  +  *£ 

*DEL 


ElZ  ;  X  0*F[1 ;  l]xj?[2  ;  1  ] )  x£[2  ;  ]  -El  1 ;  ] 

ZC,[3]^ZC’[3]  +  0=7’[1  ;  1  ]  xtf[  2  ;  1  3 

S+S+E 

SS<-SS+E*E 


SKIP  l  JUMP 
6  2  t7> 

( 3  5  p ’  *  )  ,  6  2  Tf 

JUMP:-*(NS>CNT+CNT+1  ) /AGAIN 


(I-l)p’  ’ 

>4«-  3  15  p  '  EXPONENTIAL 

OUTPUT:  (5p*  ’  )  ,  >4  [  J  ;  ] 

'QUANTITY  ESTIMATED 

(1*3) /' ACTUAL  VALUE 

' MEAN  OF  ESTIMATES  ' , 
•  0F  ESTIMATES  '  , 


EST EXTREME  VAL  ESTDIFFERENCE 


BETA  THETA  LD1 5  7D 

95  77>99’ 

’ ,  9  3  v  1  0  1.099  2.944  4.5 

95 

9  3  tS[7  ;  3+N+-NS-  ZCll  ] 

9  3  T(SS[7;]*AM)-(S[7;]*2)*tfxtf 

-  1 


. 

1  ;\ 

T' 
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[44]  'ALSO  THERE  WERE  '  %  ( vZC [ I ] ) , ’  SAMPLES  WHICH  GAVE  NO  S 

OLUTION ' 

[45]  '  ’ 

[46]  -»(  3>I+I+1  ) /OUTPUT 

[47]  ’ THE  CPU  TIME  REQUIRED  WAS  ’  , T 0 . 0 0 1 *LM I  [  2  ]  -  T 


no 


[O] 

[1] 

[2] 

■  C3] 

[4] 

[5] 
C6] 

[7] 

[8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 

[29] 

[30] 

[31] 


EVMOM \X 1 ; DEL ;XN ; P ; NS ; ND ; R ; S ; 55 ; E ; VAR ; A  V ; T ; CNT \ZC\B\ER 

R 

' THIS  PROGRAM  CALCULATES  THE  ESTIMATES  OF  BETA ,  THETA 

,  L£>7  5  ,  ’ 

1 L2795,  yUVP  L£99  [/SltfC  THE  MOMENT  ESTIMATES  FROM  THE  E 

XT RE ME' 


' VALUE  DISTRIBUTION .' 

T  7  [  2  ] 

'ENTER  THE  STARTING  POINT  XI:' 

Xl+D 


THE 


INCREMENT  SIZE  DELTA 


'ENTER 
DEL+U 
'ENTER 
ND+U 
'ENTER 
NS+  □ 

'ENTER 
URL+-  □ 

S+SS+E+5  p  ZC* 

P+1  000 Ox  +  \GEN 
R+®  1  3  19  99 

AGAIN :  XN+X 1  -  DEL*  +  /  (  ?JVZ7pl  0000  )  ° 
ZC,«-ZC,+£,flJ?<-0  =  K>r.R<-(  (  +  /XN*XN)+ND- 


THE  NUMBER  OF  DEATHS  IN  EACH  SAMPLE:' 

THE  NUMBER  OF  SAMPLES  IN  THE  SIMULATION: ' 


THE  RANDOM  NUMBER  GENERATING  SEED: 


-CNT+O 
XI , DEL 


>P 
1 )  ■ 


(  UV+(  +/XN)+ND)  *2  )  * JV 
D+ND- 1 


+ERR/ AGAIN 

B+  (Ol)  *(6x0  R)  *0.5 

S+S+E+B  ,AV+  (*£)x/?+-  0  .  57721  6+©"  l  +  *£x£>££, 
SS+SS+ExE 


-+(NS>CNT+CNT+1  )  /AGAIN 


'QUANTITY  ESTIMATED  BETA  THETA  LD1 5 

95  ‘  L£99’ 

’  ACTUAL  VALUE  '  ,  9  3  Tl  ,/? 

'MOW  <9F  ESTIMATES  ’,93 

’07?  0.F  ’,9  3  t  (  SSVNS-  1  )  -  (S*S)  *-NS*NS-  1 

OLS<7  Ttfj?/?#  WERE  ’  ,(t ZC),’  SAMPLES  WITH  NO  SOLUTION.' 
'  THE  CPU  TIME  REQUIRED  WAS  '  ,  t  0  .  0  0  1  x\jA  I  [  2  ]  -  T 


. 

. 


Ill 


[0] 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 

[29] 

[30] 

[31] 

[32] 

[33] 

[34] 

[35] 

[36] 

[37] 
[3  8] 

[39] 

[40] 

[41] 

[42] 

[43] 


FIN\P\X\K\N\NS\X2\T\NC\B\TH\NB\NTH\A\C\D\E\ DEN \F\FF\X 


FF;G;H ; ERR ; 5 ; 55 ; CNT ; R ; J 
' THIS  PROGRAM  USES  FINNEYS  METHOD  TO  ESTIMATE ' 

'BETA,  THETA ,  1575,  LD 95,  AND  LD 99.  A  FIXED ' 

’ NUMBER  OF  SAMPLES  OF  THE  SAME  SIZE  ARE  TAKEN ' 

'AND  THE  NUMBER  OF  DEATHS  IN  EACH  SAMPLE  ARE ' 

’ USED  TO  FIND  THE  MAX  LIK  SOLUTIONS . ' 

' ENTER  THE  PERCENTILES  AT  WHICH  SAMPLES  ARE  TAKEN:' 


P+ 0  .  Olx  (  i  0  )  ,□ 

'ENTER  THE  SAMPLE  SIZE:' 
K+U 


'ENTER  THE  NUMBER  OF  SAMPLES  IN  THE  SIMULATION:' 
NS+U 

'ENTER  THE  RANDOM  NUMBER  SEED:' 

DRL+-D 


X2+X*X+®P±1-P 
R+®  1  3  19  99 

S+SS+  5  p  NC+CNT+UIO+  0  x  T+UA  II  2  ] 

P+§+\  lOOOOx  (Po  .  *XK)  x  (  (  l-p)  o  .  *K- \K)  x  (  (  PP)  ,K)  p  (  \K)  IK 

Uio+i 

MA I NLOOP :  N++  /  (  (5,pX)p?(pX)pl0000)>5 

i+b+i+th+o 


AGAIN:  A<-+  /FF<-Fxl-F<-l  +  l  +  *-BxX-TH 
ERR+-  0  .  0  0  0  1  >  |  DEN+KxB*  (A*D<-+  /X*XFF)  -  ( C++  /XFF+X^FF  )  *  2 
NC+NC + ERR+E RRv  1  0  <1+1  +  1 


-+ERR  /  MAI  NLOOP 

NTH+TH+ ( +DEN) x ( ( G++ / N- K*F ) x - D- T H*C ) + ( H++ / X*N - K*F ) *C -T 

H*A 


NB+B+  ( +DEN)  xfix  (tfx/ ) -£xC 

NC+NC  +  ERR+  (NB<0  .l)v(/VB>10)v25<|  575 

-»5BB/M;1  INLOOP 

-*(  (  0  -  0 1  <  I  ( B+NB  )  -  B  )  v  (  0 . 0 1  <  I  (TH+NTH)  -TH)  )  /AGAIN 

S+S+E+B , TH+R+B 

SS+SS+E*E 

+(NS>CNT+CNT+1 ) /MAI NLOOP 


QUANTITY  ESTIMATED 
ACTUAL  VALUE 


BETA  THETA  LD1 5 


95  1599’ 

3  t  1  ,  B 


15 


A/BMA/  55  557IMvl7B5  ’,93  *5^55 

VAR  OF  ESTIMATES  ’,93  ▼ ( SS+NS- 1 ) - ( 5x5 ) +NS*NS~ 1 

t 

755  EXPECTED  NUMBER  OF  DEATHS  IS  ' ,*K*+ / l-*l+*-X 
ALSO  THERE  WERE  ’ ,(tM7),’  SAMPLES  WHICH  DID  NOT  CONV 

ERGE . ' 


THE  CPU  TIME  REQUIRED  WAS  '  ,  t 0  .  0  0  1  x[}4 1  [  2  ]  -  T 


: 

j 

s  •  s 
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[0] 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 


R<-GEN  PAR\X  1  i  DELTA  ;  NUM ;  PSUM  \R\Q\N \P\PI 
p.  THIS  FUNCTION  GENERATES  THE  PROBABILITIES  FOR 
ft  THE  FIRST  ZERO  DISTRIBUTION .  THETA  AND  BETA 
R  ARE  ASSUMED  TO  BE  0  AND  1  RESPECTIVELY 
Xl^PARll ] 

DELTA+PAR12 ] 

NUM+PARlpPAR]+ 500* ( pPAR) <3 

PSiW«-i?<-l-G<-l*l  +  *-  21 

A^l 

LOOP:  P<-QxPI+-l  -1t-1  +  *-  Ji  -DELT A  *N 
PSUM-^PSUM+P 
R<-RtP 

Q<-Qxl-Pl 

N+-N+1 

->(  ( N<NUM )  *PSUM<0  .  9  999  )  /  LOOP 


; 
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[0]  R+GENOR  PAR; XI ; DELTA ; NUM ; PSUM ; R;Q ; N ; P ; PI 

[1]  a  THIS  FUNCTION  GENERATES  THE  PROBABILITIES  FOR 

[2]  a  THE  FIRST  ZERO  DISTRIBUTION  WHEN  THE 

[3]  a  UNDERLYING  DISTRIBUTION  IS  NORMA  1(0,1) 

[4] 

[5]  DELTA<-PAR[2'] 

[6]  NUM*-PARtpPAR']  +  500x(pPAR)<Z 

[7]  PSUM+-R+1-Q+NORMAL  XI 

[8]  N+- 1 

[9]  LOOP:P*-Q*PI+l-NORMAL  XI -DELTA  x  N 

[10]  PSUM+PSUM+P 

[11]  R+R,P 

[12]  Q+Qxl-PI 

[13]  N+N+ 1 

[14]  ->(  (N<NUM)  *PSUM<0  .999  )  /  LOOP 


' 
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[03 

[  1 3 
[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 

[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 


XI  NO REV  DEL\XN\P\NS\ND\R\S\SS\E\V AR\AV \T \CNT \ZC \B \ER 

R 


' THIS  PROGRAM  CALCULATES  THE  ESTIMATES  OF  BETA ,  THETA 

,  LP75,’ 

1 LD 95,  AND  LD 99  USING  THE  MOMENT  ESTIMATES  FROM  THE  E 

XT RE ME1 

' VALUE  DISTRIBUTION  WHEN  THE  ACTUAL  DISTRIBUTION  IS  N 

ORMAL ' 

'WITH  MEAN  0  AND  VARIANCE  PI* 2/6.' 

T  J  [  2  ] 

'SPECIFICATIONS: ' 

'  XI  DELTA  SAMPLE  SIZE  SIMULATION  SIZE 

SEED' 

7372130170180  tXI  ,  DEL  ,  NUMDEA  THS  ,  SIMSIZE 

ND<-NUMDE  AT  HS 

NS+-SIMSIZE 


S+SS+-E+ 5  pZC+CNT+-0 

P+-1  0  0  0  x  +  \  GENOR  (XI,  DEL  )  *  (  6*0  .  5  )*01 
P^©  1  3  19  9  9 

j4Gj4IN:  XP*-X1  -DEL*  +  /  (  ? NDp  1  000  )  o  .  >p 
ZP^ZP  +  FFF*-0  =  yZF*-(  (  + /XN*XN)±ND-  1  )  -  (  (  AV+ (  + /  XN)±ND)  *2)  *N 


DRND-1 

-+ERR/ AGAIN 
B«-(Ol)*(6x7//?)*0 .5 

S+S+E+-B  tAV+(*B)*R+~  0 . 5  7  7  2  1  6  +  ©"  1  +  *B*DEL 


SS^SS+ExE 

-*(NS>CNT+CNT+1  ) /AGAIN 
'QUANTITY  ESTIMATED 

'ACTUAL  VALUE 
'BIAS  OF  ESTIMATES  ' ,  9 
’  IMP  OF  FFFIM^FFP  ’  ,  9 

(tZC),  ’  5^A7PIF5  //O 


5FF/1  F  HET  A  ’  LD  IS  LD 

95  LD 99’ 

?  ,  9  3  tF*-0  ,P75  ,P95  ,P99 
3  v (S+NS) -  1 ,E 

3  ?(SS+NS-1 ) - (SxS)+NS*NS-l 
SOLUTION.  CPU  TIME  '  ,  (tO.OOI* 

□i4i[2]-n 


f  f 

T  t 


* 

. 


' 


’ 

h 

' 

m 


[0]  R^ NORMAL  X  ;X  SQ  ;  P ;  REM ;  //  \ERR\  SI  G  N ;  L  ;  H 

[1]  a  THIS  FUNCTION  RETURNS  A  NORMAL ( 0 , 1 ) 

[2]  a  DISTRIBUTION  FUNCTION  VALUE 

[3]  a  ACCURACY  IS  TO  4  DECIMAL  PLACES 

[4]  2>x>~4 

[5]  //*-*<  4 

[6]  X<-X  *L*H 

[7]  SIGN+-X<  0 

[8]  XSQ+X*X*-\X 

[9]  P^(*~0.5x*S£)x(O2)*~0.5 

[10]  ERR+0  .  0001-fP 

[11]  R+-REM+X 

[12]  tf<-l 

[13]  LOOP:N<-N+2 

[14]  REM+REMxXSQ*N 

[15]  R<-R  +  REM 

[16]  ->  (  v  /  REM>  ERR )  I  LOOP 

[17]  RESIGN + ( 1- 2*SIGN) x0 . 5  +PxP 

[18]  P^(7?xL)+o  .  5x~p 
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[0] 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 
[9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 


XI  NORSIM  DEL \T \R\K\N \P\ NS \D\B\ NO  SOL ; T  H \ E ; SUM ; SSQ ; ERR 

t  I 


’ THIS  FUNCTION  CALCULATES  THE  EXPONENTIAL  ESTIMATES' 

' WHEN  THE  ACTUAL  UNDERLYING  DISTRIBUTION  IS  NORMAL 
T  J  [  2  ] 

R+®  1  3  19  99 
’ SPECIFICATIONS : ’ 

’  XI  DELTA  SAMPLE  SIZE  SIMULATION  SIZE 

SEED' 

6  3  6  2  11  0  18  0  18  0  vX 1 , DEL , NUMDEA THS , SIMSIZE , URL 

K+-NUMDE  AT  H  S 

N+SIMSIZE 


P+1000*+\GENOR(X1  ,M)f(01  )  x6*-0 . 5 
SUM+-SSQ+-SpNS^NOSOL^O 
MAI  NLOOP :  D<-DELx  0  .  5  +  +  /  (  ?£pl  000  )  °  .  >P 
B+-BEXPML  D 

NOSOL<-NOSOL  +  ERR+-(B>10  )  vB<0  .  1 
-+ERR/MAI  NLOOP 


TH<-X1+  (-f£  )  x©£x  (  “  l  +  *B*DEL)±  +  /~  1  +  *B*D 


E^B , TH+R+B 
SUM+SUM+E 
SSQ+SSQ+E* 2 
(  N>NS*-NS+  1  )  /MAI NLOOP 
' QUANTITY  ESTIMATED 

' ACTUAL  VALUE 
' MEAN  BIAS  '  ,  9 

'VAR  OF  ESTIMATES  '  ,  9 
(t NOSOL),'  SAMPLES  GAVE 


BETA  THETA  LD1 5  LD 

95  LD 99’ 

’,93  JE+-  0  ,  N1  5  ,  // 9  5  ,N  99 
3  v (SUM+N) -1 ,E 
3  t (SSQ+N- 1 ) - (SUM*2 )*N*N- 1 
M7  SOLUTION .  CPI/  7IMP  '.tO.OO 

1xD/J[2]-7’ 


t 


t 


■ 

' 
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[0]  71  SINKEV  DELTA  ;  P ;  D  \  DROP ADD  ;V 

[1]  r  GET  BASIC  PROBABILITIES 

[2]  P++\$GEN  71,  DELTA 

[  3  ]  D+-  (  ©“  l  +  *DELT A  )  +4>7 1  +  DELTA  x  1  -  n  p P 

[4]  fi  TRIM  BELOW 

[5]  DROP<-+/D<~  3 

[6]  P+DROP\P 

[7]  D+DROPi  D 

[8]  r  TRIM  ABOVE 

[9]  DROP++ /D>H .b 

[10]  D<-  (  -DROP)  i-D 

[11]  D+-~2\ ,D,  (  0 . 5x£+l4>£>)  ,  [  1  .  5]  1<J)Z?  • 

[12]  P+-  (  -  DROP)  \P 

[13]  P+-  ~2l,P,P,[1.5]  P 

[14]  r  CALCULATE  THE  TRUNCATED  PROBABILITIES 

[15]  Y3<-*-*-Z> 

[16]  Y  3+-Y  3  ^  “  1 1 Y  3 

[17]  r  FILL  OUT  THE  VECTORS  TO  4 . 5 

[18]  V<-~  3  .  1  +  0  .  2x  i  38 

[19]  ADD<-+  /DlpDl>V 

[20]  D^DyADDiV 

[21]  P+-P , ADD\ 3  8  p 1 

[22]  .  Y3+-Y3  ,  A  DD  4-  3  8  p  1 

[23]  r  PUT  INTO  DISCRETE  FORM  IE  STEPS 

[24]  Y+-D 

[25]  Y 1  <~P 

[26]  Y2<-*-*-7 

[27]  (  v  p  7 )  ,  '  ,  3,* 

*•[28]  ( (p7)  ,56)p(  (  (4xp7)  ,13)p  13  4  t7 , Y 1 , Y 2 , [ 1 . 5 ]  Y3),’,’ 


■ 

■ 


' 


[0]  XI  SINKEXP  DEL; DROP 
[  i  ]  P«-+  \  GEN  XI  ,  DEL 

[2]  D<-((*-Xl  )  +  ~l  +  *DEL)x'l  +  *DEL*~l  +  \pP 

[3]  DROP<-+/D>b 

[4]  D<-(  (-DROP)  ±D)  ,  5 

[5]  ZH-"  2  4  %D%  (  0  .  5xD+l<t>p)  4  [  l  .  5  ]  i$D 

[6]  P«-( 1-DROP) iP 

[7]  P<~ "24',P,P,[1.5]  P 

[8]  E<rl-*-D 

[9]  (tP D) , ’ ,  2, ’ 

[10]  (  (  pD)  ,  42  )  p  (  (  (  3xpD)  ,  13  )  p  13  4  tZ?,£,[1.5]  P) 


y 


■ 
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[0] 

[1] 

[2] 

[3] 

[4] 

[5] 

[6] 

[7] 

[8] 
C9] 

[10] 

[11] 

[12] 

[13] 

[14] 

[15] 

[16] 

[17] 

[18] 

[19] 

[20] 
[21] 
[22] 

[23] 

[24] 

[25] 

[26] 

[27] 

[28] 

[29] 

[30] 

[31] 

[32] 

[33] 

[34] 

[35] 

[36] 


THETA  PAR-,  XI', DEL \XI\WI\TH\ NTH \E\G\GP\ ND ; NS ; CNT ; T ; P ; S ; 

55 ; WEX ’,  N ;  EX \D \T  OT \T  0 \ A  ;  B 

' THIS  PROGRAM  CALCULATES  THE  EXACT ,  EXPONENTIAL,  AND 

EXTREME  VALUE' 

' ESTIMATES  OF  THETA  AND  THEIR  DIFFERENCE  WHEN  BETA  IS 

KNOWN  TO  BE  1 . ’ 


Xl+PARll ] 
DEL+-PAR12  ] 
ND+-PARI  3  ] 
NS+SIMSIZE 
' SPECIFICATIONS' 


'  XI  DELTA  DEATHS  SAMPLES  SEED ' 

103102100100140  tXI ,DEL, ND,NS, URL 
T+-UAI  12~] 

P-+1  0  0  0  0  *  + \ GEN  XI, DEL 
5<-55^8  p  T OT^C NT+-0 
A<~-  0  .  57721567-®”!  +  *DEL 
B+-  (  ®ND )  -  PS  I  ND 

MAI NLOOP : N^l  +  +  / (?  NDplOOOO) ° . >P 
-+(  =  /N)  / MAINLOOP 

E*-TH+X1+®{ND*~  1  +  *DEL)±+  /~  1  +  *DEL*N-  1 
WI<-+/ (  i  r /N)  o  .  <N 
WEX<-WIxEX+*-XI^Xl+DEL*l-  i  p WI 
THLOOPi  G<-  (  -  ND)  +  +  /WEX* D+-+EX  +  *T H 
GP<-(*TH)*  +  /WEX*D*D 


NT  H+-T  H  +  G^GP 

+  (25< I  NTH) /MAINLOOP 

•-*  (  0 . 0 1  <  I  (  TH+-NTH )  -  TH)  /  THLOOP 

NTH^A+Xl+DELxl -TO^ ( +/N)+ND 

T0T+-T0T  +  T0 

S*-S+E<-T  H  ,  ( TH-E )  ,E,B,  ( E-B )  ,  ( E-B  +  NTH )  ,  NTH  ,  NTH-TH 
SS+SS+ExE 

-*(NS>CNT+CNT+1 ) / MAINLOOP 


'TYPE  OF 

' ESTIMATE 

'BIAS',  9 
’ VAR  ',  9 
'CPU  TIME 


□<-0-0- 1  ’ 


EXACT  EXPONENTIAL 

EXT  VAL ' 

MLE  MLE 

MOMENT ' 

3  vS+NS 

3  v  (  SS+NS-  1  )  -  (  5><5 )  +NSXNS~  1 
' , (vO .001xCMI[2]-7) , ’  ■  OBS/DEATH' , 

T  -VNS 


UMVUE 


6  2  ?  TO 


